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FOREWORD 


'Fhe present report is one of a series of six reports, published sim 

taneously, which describe analyses and computational procedures for: 1) pre- 

diction of the in-depth response of charring ablation materials, based on one- 
dimensional thermal streamtubes of arbitrary cross-section and considering 
general surface chemical and energy balances, and 2) nonsimilar solution of 
chemically reacting laminar boundary layers, with an approximate formulation 
for unequal diffusion and thermal diffusion coefficients for all species and 
with a general approach to the thermochemical solution of mixed equilibrium- 
nonequilibrium homogeneous or heterogeneous systems. Part I serves as a 
summary report and describes a procedure for coupling the charring ablator 
and boundary layer routines. The charring ablator procedure is described in 
Part II, whereas the fluid-mechanical aspects of the boundary layer and the 
boundary- layer solution procedure are treated in Part III. The approximation 
for multicomponent transport properties and the thermochemistry model are 
described in Parts IV and V, respectively. Finally, in Part VI an analysis 
is presented for the in-depth response of charring materials taking into ac- 
count char-density buildup near the surface due to coking reactions in depth. 


The titles in the series are; 


Part I 

_ Part II 

Part III 
Part IV 

Part V 

Part VI 


Summary Report: An Analysis of the Coupled Chemically Reacting 

Boundary Layer and Charring Ablator, by R. M. Kendall, E. P. 
Bartlett, R. A. Rindal, and C. B. Moyer. 

Finite Difference Solution for the In-depth Response of Charring 
Materials Considering Surface Chemical and Energy Balances, by 
C. B. Moyer and R. A. Rindal. 

Nonsimilar Solution of the Multicomponent Laminar Boundary Layer 
by an Integral Matrix Method, by E. P. Bartlett and R. M. Kendall. 

A Unified Approximation for Mixture Transport properties for Multi- 
component Boundary-Layer Applications, by E. P. Bartlett, R. M. 
Kendall, and R. A. Rindal. 

A General Approach to the Thermochemical Solution of Mixed Equilib- 
rium-Nonequilibrium, Homogeneous or Heterogeneous Systems, by 
R. M. Kendall. 

An Approach for Characterizing Charring Ablator Response with In- 
depth Coking Reactions, by R. A. Rindal. 


This effort was conducted for the Structures and Mechanics Division of 
the Manned Spacecraft Center, National Aeronautics and Space Administration 
under Contract No. NAS9-4599 to Vidya Division of Itek Corporation with Mr. 
Donald M. Curry and Mr. George Strouhal as the NASA Technical Monitors. The 
work was initiated by the present authors while at Vidya and was completed 
by Aerotherm Corporation under subcontract to Vidya (P.0. 8471 V9002) after 
Aerotherm purchased the physical assets of the Vidya Thermodynamics Depart- 
ment. Dr. Robert M. Kendall of Aerotherm was the Program Manager and Prin- 
cipal Investigator. 
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ABSTRACT 


This report summarizes analyses and computational procedures for pre- 
dicting the transient in-depth response of charring ablation materials, either 
coupled to a nonsimilar, laminar, multicomponent, chemically- reacting boundary- 
layer computational procedure or partially decoupled through the use of con- 
vective transfer coefficients. The detailed developments are presented in 
companion documents. The computational procedure for charring ablators is 
an implicit finite difference procedure for an ablating surface material with 
several nonablating backup materials, it considers one-dimensional heat and 
mass transfer along thermal streamtubes of arbitrary cross-sectional area and 
permits a multiple-reaction model for gas decomposition and a general thermo- 
chemical surface boundary condition. The boundary-layer procedure utilizes a 
newly developed integral matrix solution procedure. It applies for general 
chemical systems, allowing rate-controlled surface reactions, and incorporates 
approximate formulations for mixture transport properties, including unequal 
diffusion and thermal diffusion coefficients for all species. Analyses are 
also presented for extending the boundary-layer computational procedure to 
include mixed equilibrium-nonequilibrium, homogeneous or heterogeneous general 
chemical systems and to include radiation absorption and emission, and for 
extending the charring ablation procedure to include char-density buildup due 
to coking reactions in depth. 
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AN ANALYSIS OF THE COUPLED CHEMICALLY REACTING 
BOUNDARY LAYER AND CHARRING ABLATOR 

SECTION 1 
INTRODUCTION 

The transient response of charring ablation materials to actual or 
simulated superorbital reentry depends upon the intimate coupling which exists 
between the internal heat and mass transfer processes, the surface phenomena, 
and the boundary layer which envelops the heat shield. In the absence of pro- 
cedures for obtaining fully coupled solutions, the conventional approach in 
the past has been to focus attention orf the in-depth charring ablator response, 
utilizing empirical correlations such as heat— of— ablation or ablation-rate— 
versus-surface-temperature relationships to provide the surface boundary con- 
dition. This method has been effective when applied to conditions which do 
not differ significantly from the test conditions at which the empirical re- 
lationships were derived, but there is no valid basis for extrapolation to 
other conditions since the highly nonlinear coupling between the various 
phenomena cannot be scaled. A somewhat more sophisticated approach has been 
to represent the boundary-layer heat and mass transfer processes by convective 
transfer coefficients while considering detailed chemical interactions and 
mass and energy balances at the surface. Because detailed surface physics can 
be retained in the formulation, this method is better suited for application 
at conditions beyond the range of available experimental data. However, it is 
still severely limited in this regard since the effects of nonsimilarities be- 
tween boundary-layer profiles cannot be treated precisely. Thus, mass addition, 
chemical reactions, and multicomponent diffusion effects can be taken into 
account only in an approximate manner, and upstream effects, thermal diffusion, 
and radiation-convection coupling cannot be considered except, possibly, 
through correlations of boundary-layer solutions. 

In the present series of reports, theoretical analyses are presented and 
computational procedures are described for predicting the one-dimensional 
transient response of charring (or noncharring) ablation materials intimately 
coupled to quasi-steady, two-dimensional, nonsimilar, laminar, chemically 
reacting boundary layers. In addition, procedures are described for obtaining 
charring ablation solutions with the boundary layer represented by convective 
transfer coefficients. The physicochemical models which are employed are 
outlined in the ensuing paragraphs, followed by an introduction to the specific 
computer codes which have been developed. 

Heat and mass transfer within the charring ablator is considered to be 
one-dimensional, but the thermal streamtubes are allowed to have arbitrary 
cross-sectional area. A general model for in-depth decomposition is considered. 
Detailed surface thermochemistry is considered, including selected rate-controlled 



reactions , and liquid-layer removal and mechanical spallation are taken into 
account through the use of a fail temperature for each candidate surface 
material. An approach for including char-density buildup due to coking re- 
actions is presented, but this has not been incorporated into the computa- 
tional procedure. 

The boundary-layer computer program applies to laminar axisymmetric 
or planar flow. No similarity approximations are imposed and surface dis- 
continuities (e.g. , due to change of ablation materials) are allowed. The 
procedure applies to any chemical system, considering equilibrium with the 
exception that selected species can be considered as frozen in the boundary 
layer while undergoing rate-controlled reactions at the surface. The bound- 

l^ycr procedure also considers unequal diffusion and thermal diffusion 
coefficients for all species through the use of convenient approximations to 
these coefficients. Additional theoretical developments which have been 
made but which have not been incorporated into the computer program include 
a general mixed equilibrium-nonequilibrium, homogeneous or heterogeneous 
chemical procedure and a model for radiation absorption and emission in the 
boundary layer. 

The convective transfer coefficient approach utilizes the same charring 
ablation computational procedure (with its general surface thermochemistry, 
decomposition model, and thermal streamtube approach) and is less restrictive 
in one sense in that the boundary layer can be turbulent as well as laminar. 
With regard to other boundary- layer phenomena, it permits consideration of 
unequal diffusion coefficients, but neglects thermal diffusion and includes 
nonsimilar effects only if they have been determined from correlations of 
boundary-layer solutions. 

The boundary-layer computational procedure utilizes an entirely new 
numerical solution procedure which was developed specifically for the present 
problem. It has come to be known as an integral-matrix method because of the 
way the problem is formulated and solved. The charring ablation and chemistry 
routines and the convective transfer coefficient approach are extensions of 
procedures which have been under continued development by the present authors 
during the past several years. 

The following computer programs have been developed and are described 
in this series of documents: 

1. Charring material ablation (CMA) program 

2. Boundary-layer integral-matrix procedure (BLIMP) program 

3. Surface . thermochemistry programs based on convective transfer 
coefficients: Aerotherm chemical equilibrium (ACE) and Aerotherm 
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chemical equilibrium with selected surface reaction kinetics 
(ACE-KINET) 

4. Coupled ablator/boundary layer/environment (CABLE) program 

The CMA program can be operated independently for obtaining the in-depth 
response of charring materials for assigned ablation rates and surface tem- 
peratures. The ACE and ACE-KINET programs can be operated independently 
to compute steady-state ablation of arbitrary material-environment combina- 
tions or can be used to provide tables of surface-thermochemistry information 
as input to the CMA program for transient charring ablation calculations, 
again for arbitrary material-environment combinations. The BLIMP program 
can be operated independently to provide boundary- layer solutions for a va- 
riety of coupled, partially coupled, or uncoupled steady-state ablation sur- 
face boundary conditions. Finally, the CABLE program calls upon the BLIMP 
and CMA programs as subroutines to provide fully coupled transient charring 
ablation and boundary-layer solutions . 

The analyses and computational procedures for characterizing the laminar 
chemically-reacting boundary layer, for predicting the in-depth transient 
response of charring materials, and for evaluating the chemical state are 
summarized in Sections 2 through 4 and are presented in detail in Parts III, 
II, and V, respectively. The approximations for multicomponent transport 
properties are summarized in Section 2.1.2 and reported in detail in Part 
IV, the approach for characterizing charring ablator response with m-depth 
coking reactions is summarized in Section 3.5 and described in detail in 
Part VI, and the radiation model is summarized in Section 2.1.3 and presented 
in an Appendix to Part III. The coupled boundary- layer and charnng-ablation 
procedures (i.e., CABLE and CMA/ACE and CMA/ACE-KINET combinations) are des- 
cribed in Section 5. Summary, conclusions, and recommendations are pre- 
sented in Section 6. 


SECTION 2 

ANALYSIS AND COMPUTATIONAL PROCEDURES FOR LAMINAR 
CHEMICALLY- REACTING BOUNDARY LAYERS 

The boundary layer which envelops an ablating heat shield during super- 
orbital reentry is intimately coupled with the transient ablation processes. 
In addition: 

1. It may be laminar, transitional, or turbulent on different parts 
of the body and at various flight conditions. 
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It may be highly nonsimilar, especially if there are changes 
of ablation materials. 

3. The surface material may react chemically with the air (or other 
planetary gas) , change phase, and/or be removed mechanically by 
spallation or liquid-layer runoff. 

4. Chemical reactions will generally also occur throughout the 
boundary layer. 

5. The homogeneous and heterogeneous reactions may be kinetically 
controlled . 

6. Incident radiant energy may be absorbed at some wave lengths 
and emitted at other wave lengths . 

7. The molecular species in the boundary layer will be governed 
by different diffusion coefficients relative to the other 
molecular species which are present. 

8. Thermal diffusion can be significant, especially if there are 

wide variations of molecular weight (e.g., when hydrogen is present 
in the boundary layer as a result of the decomposition of a charring 
ablation material) . 

9. An entropy layer may be present. 

10. At very high entry velocities, the inviscid flow field may be 
nonadiabatic as a result of radiation cooling. 

In the present report, only the laminar boundary layer is considered in 
detail. Otherwise, all of the features listed above are considered in the 
theoretical analysis. However, Items 5, 6 and 10 are not presently allowed in 
the computational procedure, with the exception (with regard to Item 5) that 
selected species can be considered to be frozen in the boundary layer and to 
react with finite rates at the surface. The major features contained in the 
theoretical analyses are summarized in Section 2.1. The numerical solution 
procedure is introduced in Section 2.2. The computer program, designated 
BLIMP for boundary layer integral matrix procedure, is described briefly in 
Section 2.3. The analysis and numerical procedure are presented in consid- 
erably more detail in Part III of this series, while a user's guide to the 
BLIMP program is presented in Ref. 1. 

2.1 THEORETICAL ANALYSIS 

The boundary— layer analysis which has been developed is based on the 
following model: 
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1. Laminar axisymmetric or planar flow with all nonsimilar terms 
retained and discontinuous wall conditions allowed (e.g., change 
of heat-shield material) 

2. Multicomponent transport properties (including unequal diffusion 
and thermal diffusion coefficients for all species, based on newly 
developed approximations for these coefficients) 

3. Coupled radiation absorption and emission using a conventional one- 
dimensional model 

4. General boundary conditions, including mass and energy balances at 
the wall and the presence of an entropy layer 

5. Mixed equilibrium- nonequilibrium, homogeneous or heterogeneous 
chemistry for general chemical systems. 

The first four features are discussed in Sections 2.1.1 through 2.1.4, respec- 
tively.. The chemical model is discussed in Section 4. 

2.1.1 Conservation Equations for Nonsimilar, Laminar, Plan ar or 
Axisymmetric Boundary Layers 

Similarity approximations have been often applied in order to 
help simplify complex boundary-layer problems. This involves a transforma- 
tion to a new coordinate system (T],§) from the original (y,s) system where 
y and s are the normal and streamwise coordinate, respectively. The simi- 
larity transformation is successful if the ^-variations of functions of the 
dependent variables vanish or become of negligible importance, in which case 
the partial differential equations become ordinary differential equations. 

The most popular transformation, known among other names as the Levy-Lees 
similarity transformation (Ref. 2) , is given by 

pdy (i) 


f = f up M r 2K ds (2) 

b J e H e H e o 

o 

where the subscript e refers to the boundary- layer edge, p is the density, 
u is the streamwise velocity, pi is the viscosity, r Q is the local radius 
of the body measured normal to the body centerline, and k is zero for 
planar bodies and unity for axisymmetric bodies . 


T] = 


r H u 
o e 

( 2 ?) h 


J 
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It can be shown that the similarity transformation is valid at the stag- 
nation point since the terms involving ^-derivatives vanish at the stagnation 
point (§ = 0) , and is valid within certain restrictions on the surfaces of 
flat plates, wedges, and cones. These restrictions include; 

1* No streamwise variation of surface properties 

2. Streamwise variation of boundary- layer-edge properties in accordance 
with specific relations 

3 • An inverse-square-root type variation with ? of surface mass- 
transfer 

4. Transport properties functions of tj only 

5. Edge conditions functions of ? only (i.e., no entropy-layer 
effects) 

Furthermore, the similarity assumption is not valid on bodies of arbitrary 
shape . 

The importance of nonsimilar effects is illustrated in Fig. 1 where 
constant blowing into an incompressible boundary layer over a flat plate is 
terminated two feet from the leading edge. (This solution was generated with 
the BLIMP procedure to be described later in this report.) In this example, 
only Item 3 in the above listing is violated. All results shown in Fig. 1 
are for a nonsimilar boundary layer. The point is that a similar solution 
downstream of injection would show immediate recovery of the wall shear to 
the asymptotic value, whereas the wall shear is still 13 percent below this 
value two feet after blowing is terminated. 

It is apparent from the above discussion that a similarity assumption 
is overly restrictive for a general solution procedure. Therefore, in the 
present analysis, although a similarity transformation is employed because 
of the normalizing benefits therefrom derived, all nonsimilar terms are re- 
tained. That is, no similarity assumptions are made. 

The transformation which is employed is similar to the Levy-Lees trans- 
formation (Eqs . (1) and (2)) , with two exceptions; 

1. A coordinate stretching parameter is introduced, 

n = Tj/a H (§) , T =* ? (3) 

where a R (?) is a dependent variable determined during the course 
of the calculation. The purpose of this supplemental transformation 
is to constrain the maximum and minimum values which the transverse 
coordinate can assume in a nonsimilar boundary layer. 
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2. The subscript e refers to the wall streamline value for the 
inviscid solution, which differs from the boundary-layer-edge 
condition in the presence of an entropy layer. 

Application of this transformation yields the following boundary— layer 
conservation equations. 


Momentum equation 


ff 1 


rcf-i 
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Energy equation 

fH i + [ -q a + <$] 
Elemental species equations 
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where the prime denotes partial differentiation with respect to rj, f is 
the stream function defined by 



o 


such that f‘ = a„u/u , f is given by 

rt e W 


f 

P w v w tota ^- mass flux into the boundary layer, is the total 

enthalpy, is the mass fraction of base species k independent of 

molecular configuration, C = p|i/pjj , g = 2d in u /d in T, <f>. is the 
source term for base species k which arises in nonequilibrium systems 
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(see Section 4) , q is the diffusive heat flux, q is the net one dimen 
sional radiant heat flux, is the diffusive mass flux of base species k, 

and the asterisk signifies normalization by division by a*, defined by 

a * s P e u e^e r o /(2l) ^ (9) 


Expressions for q a and including approximate formulations for unequal 
diffusion and thermal diffusion coefficients are presented in Section 2.1.2, 
whereas a relation for q^ based on a one— dimensional model is presented in 
Section 2.1.3. 


2.1.2 Diffusive Fluxes in a Multicomponent Boundary Layer Based on 
Approximations for Unequal Diffusion and Thermal Diffusion 
Coefficients 

Consideration of unequal diffusion and thermal diffusion coefficients for 
all species adds considerable complexity to the boundary-layer solution if 
these coefficients are treated in a precise manner (i.e., first-order kinetic 
theory for the binary diffusion coefficients and second-order kinetic theory 
for multicomponent thermal diffusion coefficients) . On the other hand, the 
assumption of equal diffusion coefficients for all species and the neglect of 
thermal diffusion (or the use of constant thermal diffusion factors) are 
overly restrictive. Therefore, convenient approximations for these coef- 
ficients have been developed and introduced into the boundary-layer equations. 

The approximations per se are described in detail in Part IV and are briefly 
reviewed in this section. Expressions for diffusive heat and mass fluxes in- 
corporating these approximations are developed in Part III and summarized herein. 

The approximation for binary diffusion coefficients* is of the form 




D 

F i F j 


( 10 ) 


where D is a reference diffusion coefficient and F^ is a diffusion factor 
for each species in the mixture. Thus, the v (v — 1.) /2 diffusion coefficients 
pertinent to a molecular set of v species are replaced by D and v dif- 
fusion factors F^.**The primary advantage of this approximation is that it 
enables explicit expression of the diffusional mass flux of species i, j^, 
in terms of gradients and properties of species i and of the system as a 
whole. This, in turn, permits the Shvab-Zeldovich transformation to be made 
without introducing the concentration-dependent multicomponent diffusion 
coefficients, j • This transformation reduces the number of species 

*Such a bifurcation approximation was first proposed by Bird (Ref. 3) . 

**Thus, the procedure is exact for binary or ternary systems and amounts to a 
correlation of diffusion coefficient data for larger systems. 
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conservation equations from one for each of the molecular species (e.g., 20 
to 70) to one for each element (e.g., 3 to 5), A thorough study of the accu- 
racy of the approximation, presented in Part IV, shows that diffusion coef- 
ficients obtained using this technique are generally within a few percent of 
the values predicted by kinetic theory. 

It is convenient to consider to be unity for some reference species 

such as molecular oxygen. Then the D is not too unlike the self-diffusion 
coefficient for that species. The for the other species in the system 

of interest are then obtained by means of a least-squares fit of actual dif- 
fusion coefficient data. The might be expected to depend upon tempera- 

ture, pressure, and concentrations. However, the F^ are, in fact, inde- 
pendent of pressure, since the pressure dependence of can be fully ab- 

sorbed into the D. Furthermore, satisfactory accuracy can be obtained while 
considering the F^ to be independent of temperature and concentrations. 
Thus, the F^ can be considered as a set of constants for a given set of 
molecular species. in addition, it has been found that the F^ for a given 
species i will in general be nearly the same when species i is considered 
as part of two different sets of molecular species, being given approximately 
by = (77^/26) 0 * 4 6 1 where 77^ is the molecular weight of species i. 

The approximation for multicomponent thermal diffusion coefficients, 

T 

, is based on a correlation of binary thermal diffusion data and a gen- 
eralization to multicomponent systems, satisfying the requirement that 

Z T T 

D i = recognizing that is independent of fluxes, and assuming 

that thermal diffusion of species i behaves as though it were in a system 
of species i and a species representative of the mixture as a whole. The 
resulting equation is given by 


c.pDfj 

D i - i hr <z i - K i> 


where is mass fraction of species i, z^ is defined as 


in k . 

Z. = = 7 “ 

1 ^2 F i 


p, and |i ? are system properties defined by 
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i=l 
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( 12 ) 


(13) 
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with x. the mole fraction of species i, and c, an empirically deter- 

i ^ 

mined parameter which was found from the correlation of binary thermal dif- 
fusion data to be nearly constant (-0.5). This approximation permits con- 
sideration of multicomponent thermal diffusion coefficients for all molecular 
species with very little computational effort. 

The approximation for binary diffusion coefficients is also introduced 
into approximate expressions for multicomponent viscosity and thermal con- 
ductivity of the Sutherland-Wassil jewa type (see Ref. 4) to simplify them 
still further. For example, the resulting expression for viscosity is given 
by 

v 

“ ■ 

i=l 

where c ^ and C 2 can be considered as constants. Since the F ^ and D 
can be determined in advance for a particular chemical system by a correla- 
tion of diffusion data from any desired source, it follows that reasonably 
accurate descriptions for the mixture transport properties can be utilized 
without introducing any molecular-interaction models into the input format 
of the boundary- layer program. The development of this and the other approxi- 
mate expressions for multicomponent transport properties is presented in 
Part IV. 

The diffusive heat and elemental mass fluxes for a multicomponent bound- 
ary layer are substantially simplified by introduction of the approximations 
for binary diffusion coefficients and multicomponent thermal diffusion coef- 
ficients. The diffusive mass flux of species i 

Z! + {Z ± - K ± ) 


n « U - 


a u Sc 


is given by 




(15) 


x iV c i F i 

C 2 X i F i + ^l 


(14) 


with Sc a system property defined by 




Sc = — 

pDp 2 

The diffusive flux of the k th base species is similarly given by 


( 16 ) 
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where C 

P 


is the frozen specific heat defined by 


C 

P 



K.C 
1 p i 


(19) 


with C the specific heat of species i, Pr is the Prandtl number 

j^i 

defined by 



( 20 ) 


with X the thermal conductivity of the mixture, h is the static enthalpy, 
h. is the static enthalpy of species i, and K, C , and U, are defined 

X P J 

by 


V 

r — * 


V 
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( 21 ) 


In the absence of thermal diffusion and with the diffusion coefficients 
assumed equal, the can be set equal to unity. Then the = K., 

ft = h, C p = Cp, and the expression for reduces to Fick's law. 
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Finally, it is shown in Part IV that the last term on the right-hand- 
side of Eq. (15) is small compared to the other terms in the equation. Thus, 
with fair approximation 





_ c 


( 22 ) 


and the become the driving potentials for diffusion. This observation 

serves as the basis for the development of convective transfer coefficients 
including unequal diffusion coefficients for all species. These relations 
are also presented in Part IV. However, the more precise Equation (15) is 
utili-'cd in the boundary- layer computational procedure. 

2.1.3 Radiant Heat Flux in an Absorbing and Emitting Boundary Layer 

As indicated previously, the boundary-layer energy equation contains a 
term representing the net radiative flux which must be known at every point at 
which the boundary-layer equations are evaluated. This term is important at 
lunar reentry velocities (at maximum heating conditions it is comparable to the 
diffusive flux), and it becomes dominant at still higher reentry velocities. 

A model for the radiant heat flux is presented in Appendix E of Part III and 
is summarized in this section. This model is not presently included in the 
computational procedure, but is seen as the eventual replacement for the pre- 
sent radiation model which assumes that the incident radiation passes, unatten- 
uated, through the boundary layer. 

The radiation flux term accounts for the net energy extracted from the 
local radiation field and added to the internal energy of the boundary-layer 
gases. The energy in the radiation field originates, for the most part, from 
emission in the hot shock— layer gases and subsequent transmission into the 
boundary-layer region. Scattering need not be considered since it does not con- 
tribute measureably unless particulate matter is present. The atoms, ions, 
and molecules present in the boundary layer absorb energy from the radiation 
field and distribute it among their internal energy levels. This is a highly 
frequency— dependent process, and any meaningful analysis must allow for it. 
Finally, the excess energy in internal energy levels is redistributed by inter- 
particle collisions. This is an efficient process, and when local thermo- 
dynamic equilibrium exists (as is assumed here), it occurs instantaneously. 

In forming the flux integrals, it is consistent and adequate to employ the 
conventional one-dimensional approximation, viz., local net radiant heat flux 
a function of y only. The resulting expression for the radiant heat flux is 
similar to that of Cess (Ref. 5) except for a change in the edge boundary 
condition. The present model considers a boundary layer of finite thickness 
with an angular-dependent incident flux; whereas, Cess considered an absorbing 
media which extends to infinity. The resulting expression for the net radiant 
heat flux toward the surface at the nodal point i in the boundary layer is 
given by 
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Where v is the frequency and where the spectral flux, q is given by 
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In Eq. (24), is the optical depth defined by 

i 

Y, 


J, ~ f PH v dy 
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(25) 


P* v is the mixture absorption coefficient, (T) is the Planck function 

defined as 


D /rp\ s Zh V 3 

Bv ,2 , hv/kT 


c 3 (e v/ -1) 


(26) 


with the symbols h,c and k having their usual definitions, E n (t) is the 
exponential integral defined by 


E n (t) 3 J |j n_3 e t/M du 
0 


(27) 


e is the hemispheric surface emissivity, and 


,2n 
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0 < 0 tt/2 


(28) 


14 



with I ( 9 # cp) the angular-dependent specific intensity at the edge of the 

boundary 6 layer , 9 the angle between the normal to the surface and the direction 
of incident radiation, and cp the angle between a reference line in the 
surface and the trace of the incident flux on the surface. 

The mixture absorption coefficient, p* v , is obtained by summing (at 
a fixed frequency) the contributions from all atomic, ionic and molecular 
species in the boundary layer. Thus, given by 



j 


The calculation of q^ proceeds as follows. A set of frequency points 
is selected so that the continuum spectra is accurately represented. Thus, 

20 to 30 points are usually adequate to represent the 0.25 ^ hv £ 20 ev . 
frequency range of interest. Using the given concentration and temperature 
distributions, the mixture absorption coefficients and the Planck function are 
calculated at each spatial nodal point for each of the frequencies. The in- 
tegrals required for each of the optical depths (Eq. (25)) and the spectral 
fluxes (Eq. (24)) are obtained. This" continuum flux" is obtained by inte- 
grating the spectral flux over frequency (according to Eq. (23)). 

Atomic and ionic lines are usually important in reentry problems, in 
which case, a correction to this continuum flux is required. The calculation 
of the line contribution to the total flux proceeds as follows. The frequency 
range of interest is divided into increments, each of which includes a group 
of important lines. Each frequency increment is small enough so that the 
frequency variation of the continuum spectra can be neglected. Other than 
this requirement, the frequency increments selected for the lines are unre- 

to the frequency points used for the continuum flux calculation; further- 
more, the frequency increments for the lines do not have to be connected or 
have to span the entire frequency range. A set of frequency points is selected 
for each line group so that the rapid variation of the line spectra is accurately 
described. This usually requires about 20 to 30 points per line group. The 
contribution from each line group is calculated using the same equations and 
procedures described in the preceding paragraph and then combined with the 
continuum flux to obtain the total flux. 

2.1.4 Generalized Boundary Conditions 

The wall boundary conditions of interest admit the addition of chemically 
active species arising from the pyrolysis of an internally decomposing mate- 
rial, surface combustion or phase change, and liquid-layer removal. In this 
case, the surface mass flux, surface enthalpy, and surface elemental mass 
fractions are supplied through surface chemistry considerations and energy 
and elemental mass balances. 
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The surface energy and elemental mass balances can be supplied by tran- 
sient internal conduction solutions such as those described in Section 3. 

The procedure for accomplishing this will be discussed in Section 5. The 
resultant eguation for the surface energy balance is given by 


m*h + 

g g 


m*h 
c c 
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m* h . - 
r i 1 


(p v ) *h 
w w' w 
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- q; 


cond 


= 0 (30) 


where the asterisk signifies normalization by division by a* (Eq. (9) ) , m 
is the mass flow rate per unit area and hg the enthalpy of gas which enters 
the boundary layer without phase change at the surface (e.g., pyrolysis gases), 
m c is the mass removal rate per unit area and h c the enthalpy of surface 
material (e.g., char) by chemical reactions or phase change, m r is the mass 
removal rate per unit area and h^ the enthalpy of surface material (e.g., 
char) in the condensed phase (e.g., by melting with subsequent liquid runoff 
or by mechanical spallation) , h^ is the enthalpy of the gas phase at the 
wall, q*^ is the normalized diffusive heat flux away from the wall (Eq. (18) 
evaluated at the wall), q r ^ is the net radiative flux to the wall (Eq. (23) 
evaluated at the wall), and q cond = X w (dT/dy)^ is conduction into the sur- 
face material (with X^ the thermal conductivity of the surface material) 

The elemental mass balances are given by 
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( 31 ) 


where the subscripts g, c, r and w and the asterisks have the same meaning 
and is the normalized diffusive net mass flux of elemental species k 

away from the wall, given by Eq. (17) evaluated at the wall. 

The boundary- layer edge condition is allowed to be nonisentropic such 
as results from an entropy layer or radiation cooling of the inviscid flow 
field. The former is accomplished by defining the reference condition as the 
f = 0 streamline of the inviscid flow field and specifying the edge boundary 
condition as a function of f as well as § and time. This is discussed 
in Part III. 


2.2 NUMERICAL SOLUTION PROCEDURE 

A boundary-layer solution procedure has been developed for the problem 
described in Section 2.1, that of a nonsimilar, laminar, multicomponent, 
boundary layer, with general equilibrium or nonequilibrium chemistry, radia- 
hion absorption and emission, and a transient charring ablation wall boundary 
condition . 
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As discussed elsewhere in this report, these features are not all in 
eluded in the currently operational computer code; however, it is signifi- 
cant that the numerical procedure has been formulated with this ultimate 
problem in mind and that features presently not included can be added in a 
convenient manner. The solution procedure, which has come to be known as 
an integral-matrix method, is described in detail in Part III and is 
discussed summarily in this section. 

In developing the laminar boundary- layer solution procedure, the empha- 
sis was on achieving a procedure adaptable to this general environment. It 
was also required that the procedure have the versatility to treat a variety 
of boundary conditions. Additional goals of the procedure were simplicity, 
accuracy, and computational speed. Simplicity relates to problem formulation 
and is probably best measured (negatively) by the amount of judgement re- 
quired during the solution procedure and by the amount of algebra required to 
achieve formulation. Accuracy is meant to imply that the procedure will have 
the capability, in a practical limit, of yielding exact solutions. 

In developing this procedure it was desirable to consider the characteris- 
tics of existing techniques. The iterative initial value approach seemed 
inappropriate since the inclusion of radiation with its integro-di f ferential 
character can not be accomplished in a particularly convenient fashion. In 
addition. Refs. 6 and 7 indicate a rather complex convergence process. The 
11 accuracy" requirement eliminates several methods such as simple integral 
methods, perturbation solutions and semi-analytical methods. Explicit finite 
difference methods tend to require excessive computational time. Of the re- 
maining solution procedures, three types may be considered; implicit finite 
difference (e.g., Ref. 8), matrix (e.g., Ref. 9), and integral relations 
(e.g.. Refs. 10 and 11). 

In light of the goals set for the procedure to be adopted, certain spe- 
cific requirements seemed appropriate. In particular, minimization of the 
number of entries into the conservation equations required to obtain a solution 
was judged to be of prime importance as a consequence of the relatively large 
times associated with state calculations for a general chemical environment. 

In the streamwise direction, large steps are necessitated by the desire to 
couple the boundary-layer procedure to a transient internal-conduction or 
ablation solution. 

For a given accuracy, the number of entries into the conservation equations 
necessary for solution in the surface normal direction is controlled primarily 
by the nature of the functions which relate the dependent variables (and their 
derivatives) to the independent variable. Thus the continuous functions typi- 
cally used in integral relations approaches require fewer entries than the 
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discontinuous functions implied by most finite difference approximations. in 
order to permit relatively flexible profiles, sets of connected cubics were 
selected to represent enthalpy, velocity, and concentration parameters. The 
first and second derivatives of these cubics were made continuous at the con- 
necting points. These spline functions, as they are commonly known, are 
conveniently supplied through truncated Taylor series expansions for f, H 
K k and their derivatives in terms of their higher derivatives at the same 
and the neighboring nodal points. 

If the general integral relations approach is followed, weighting func- 
tions must be selected. in the present study this selection was based pri- 
marily on the complexity of the resultant algebra. Based on studies (dis- 
cussed in Part III) using Dirac delta weighting functions (i.e., a differen- 
tial approach* ) and step weighting functions similar to those used by Pallone 
(Ref. 10) indicated , when other aspects of the procedure were unchanged, no 
apparent superiority with regard to accuracy or stability. Because all of the 
complexities introduced by the generalization of the thermodynamic and 
transport properties of the system occur within divergence terms, square-wave 
weighting functions produce markedly simpler algebra and, consequently, were 
adopted for the present procedure. 

In the past when relatively large spacing in the streamwise direction 
has been desired, iterative procedures have generally been used to assure 
accuracy and stability. In many instances (e.g.. Refs. 7 and 9) these pro- 
cedures have treated the solution in a manner resembling those used for simi- 
lar solutions but with the addition of finite difference representations for 
the nonsimilar terms, a procedure which eliminates the necessity of special 
star ^ing techniques. Using this basic approach, the specific treatment adopted 
in the current method follows most closely the matrix procedure used by Leigh 
(Ref. 9) wherein the iteration is a consequence of the solution of a set of 
linear and nonlinear algebraic relations . Whereas a special successive approx- 
imation procedure was used by Leigh, the general Newton- Raphs on technique is 
used in the present procedure. This technique results in linearized coupling 
between all relations required to characterize the boundary layer, and thus 
assures a more rapid and stable iterative convergence. In addition, coupling 
to a transient conduction solution becomes straightforward, and features such 
as nonequilibrium chemistry and gaseous radiation can be conveniently added. 


*This correspondence is pointed out by Dorodnitsyn (Ref. 11) . 
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An efficient method for solution of this matrix equation is utilized which 
takes advantage of the linear Taylor series expansions and of the zeros in 
the matrix. This will he discussed in Section 2.3.1. 

Several computer results are presented in Part III which demonstrate 
that the integral matrix method is capable of yielding accurate numerical 
solutions with a minimum number of entries into the conservation equations 
(7 to 11 spline segments have consistently yielded 3 to 4 place accuracy) . It 
is flexible and versatile as the number and spacing of spline points can be 
varied at will and the number and type of chemical elements and molecular 
species can be selected arbitrarily. Finally, it can be applied to radiation 
problems since it considers simultaneously all points across the boundary 
layer at a given streamwise station. 


2.3 BOUNDARY LAYER INTEGRAL MATRIX PROCEDURE (BLIMP) COMPUTER PROGRAM 

The BLIMP program is a computer code for solving a nonsimilar laminar 
boundary layer utilizing the integral matrix solution procedure introduced 
in Section 2.2 and described in detail in Part III. The procedure also serves 
as a subroutine for coupling the boundary layer to a transient charring abla- 
tion routine (to be discussed in Section 5) . Thermal diffusion and unequal 
diffusion are included using the approximations introduced in Section 2.1.2 and 
described in Part IV. The boundary layer can be considered as nonreacting, as 
equilibrium, or as a mixed equilibrium-frozen system with specific rate- 
controlled heterogeneous (surface) reactions or surface catalyzed reactions. 

The nonreacting boundary-layer option utilizes rather general laws for speci- 
fication of enthalpy, viscosity, and Prandtl numbers as functions of temperature. 
The other options utilize the chemical procedures introduced in Section 4 and 
described in detail in Part V. (The generalized nonequilibrium procedure intro- 
duced in Section 4 and described in detail in Part V and the radiation model 
introduced in Section 2.1.4 and discussed in detail in Appendix E of Part III 
are not presently included.) The computational procedure is discussed briefly 
in Section 2.3.1, the most significant operational considerations are sum- 
marized in Section 2.3.2, and a sample problem solution is presented in Section 


2.3.3. 

2.3.1 computational Procedure 


The integral-matrix solution procedure outlined in Section 2.2 is uti- 
lized in the BLIMP program. In essence, this involves simultaneous solution 
at each streamwise station of (7 + 3 K)N + 1 linear and nonlinear algebraic 
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is the boundary-layer thickness normalizing parameter, f is the stream func- 
tion, H t is the total enthalpy, is the mass fraction of element k, 

the subscript i refers to the spline point, and the prime refers to partial 
differentiation with respect to r\, the normal boundary-layer coordinate.* 

The boundary-layer conservation equations, integrated with a weighting factor 
of unity between neighboring points and zero elsewhere, provide (2 +K) (N - 1) 
of the equations, the Taylor series expansions which express, for example, 

f i+l in terms of f i , f i' q 1 , f|", and f«; 1# provide (5 + 2K) (N - 1) 
equations, and the boundary conditions and constraint provide the remain- 

ing 8 + 3K equations. 

The solution of these equations is achieved by use of the general Newton- 
Raphson procedure. Thus, the nonlinear equations (the boundary-layer conser- 
vation equations and some of the boundary conditions) are linearized with re- 
spect to the primary dependent variables, and the errors introduced by the 
linearization are driven toward zero by iteration. This yields a matrix of 
equations of the form 



j 


d E 
n 

av. 

j 


AV. 

J 


- ERROR (E ) 
' n' 


(32) 


where E n represents the n th equation, signifies the j th primary de- 

pendent variable (f^, q, ...), ERROR (E ) represents the error in the n th 
equation resulting from the previous iteration, and AV- is the correction 
to be added to the variable Vj for the next iteration. in order to reduce 
the nonlinear equations to the form of Eqs . (32) , it is necessary to express 
the corrections for all other dependent variables in terms of corrections 
on the primary variables. For example, the correction on density is expressed 
as 


Ap 


i 




(33) 


since the pressure is constant across the boundary layer. The derivatives 
with respect to h and K k are state properties determined at each rp by the 
chemistry routine. The Ah^ is given by 


*The streamwise derivatives do not appear as primary variables as they are 
expressed in terms of local conditions and known upstream conditions by 
means of two- or three-point finite-difference relations. 
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(34) 


Ah. 


= AH_ “ 
T i 


a 3 i 
e 


Af ! 
1 

f ! 

1 


Aa 


H 


The corrections on primed quantities (e.g., Ap|) are considerably more com 
plicated and introduce second derivatives. 

Direct inversion of the entire matrix of (7 + 3K) N + 1 equations during 
each iteration would rapidly become unwieldy as the number of spline points 
and/or elements is increased. Therefore, it is significant that the solu- 
tion procedure takes advantage of the fact that the majority of the equations 
are inherently linear and that the matrix is sparse in an ordered manner. 

The first step in this process is to utilize the Taylor series expan- 
sions and linear boundary conditions to express some of the corrections 
(termed "linear corrections") in terms of the remaining corrections (termed 
"nonlinear corrections") . These are then substituted, in effect, into the 
nonlinear equations, thereby reducing the matrix to (K + 2) (N — 1) +3 equa- 
tions in terms of (K + 2) N + 3 "nonlinear corrections." (A total of K + 2 
nonlinear wall boundary conditions are not included at this point since it 
is more convenient for some coupled problems to introduce them after the 
major matrix inversion.) This matrix equation is then reduced further in 
terms of K + 2 "reduced nonlinear corrections " (Af w » AH T W and the AK k w ^ 
by inversion of a (K + 2) (N - 1) + 3 matrix. The Af w , AH T and 

AK n are then determined from the K + 2 remaining wall boundary 
-- w 

conditions . 

The sequence of events in the actual calculational procedure is summa- 
rized in Fig. 2. This illustrates the computations which are performed for 
every new case, for each new time, and for each new station. For each sta- 
tion, there is a master iteration as indicated in the figure. For each itera- 
tion, the calculation proceeds through the boundary layer (for the given 
streamwise station) from the wall to the boundary— layer edge. This is fol- 
lowed by the major matrix inversion, which expresses the "nonlinear correc- 
tions" in terms of the "reduced nonlinear corrections." The reduced non- 
linear corrections are then provided by surface considerations, as mentioned 
previously . The "nonlinear corrections" and finally the "linear corrections 
are then evaluated, and the primary variables are corrected. The iteration 
is completed when the corrections are sufficiently small that the errors in 
all linear and nonlinear equations are acceptably small. The results for 
the streamwise location are then printed out (including all desired nodal 
data) and the solution proceeds to the next time, or to the next station (in 
which case it returns to the first time) , or to the next case. 
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Figure 2. Schematic of BLIMP program. 
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2.3.2 Operational Considerations 


In this section, such operational considerations as program options, 
input requirements and output capabilities, storage requirements, and compu- 
tational times are discussed briefly. 

2. 3. 2.1 Program options 

The BLIMP program has been designed for versatility and generality with- 
in the constraint of the physicochemical model which has been adopted. There- 
fore, the user has many options at his disposal. The ablation material, 
environmental gas, molecular species, and number and spacing of spline 
segments can be chosen arbitrarily. In addition, the more significant 
options are as follows; 

1. Body shape: 

(a) axisymmetric blunt (e.g., sphere cone), 

(b) axisymmetric sharp (e.g., cone), 

(c) planar blunt (e.g., leading edge), 

(d) planar sharp (e.g., wedge). 

2. Treatment of upstream effects: 

(a) nonsimilar boundary layer with two- or three-pomt difference 
representations of upstream information, with possible dis- 
continuities , 

(b) similar boundary layer. 

3. Chemical model: 

(a) nonreacting (homogeneous) boundary layer with variable 
properties , 

(b) equilibrium boundary layer, 

(c) mixed equilibrium- frozen boundary layer with rate-controlled 
surface reactions or surface catalyzed reactions , 

4. Transport properties; 

(a) thermal diffusion and unequal diffusion coefficients, 

(b) unequal diffusion coefficients but neglecting thermal 
diffusion, 

(c) equal diffusion and neglecting thermal diffusion. 

5. Surface boundary condition: 

(a) specified wall enthalpy (or temperature) , wall total mass 
flux (or f ) and elemental concentrations, 

(b) specified wall component mass fluxes (i.e., char, pyrolysis gas 
and edge gas) and wall enthalpy (or temperature) , 
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(c) specified wall component mass fluxes with surface equilibrium, 

(d) coupled steady— state wall mass balances and either wall energy 
balance or assigned wall temperature (providing component 
fluxes and surface equilibrium) . 

(e) coupled mass and energy balance at the wall as provided by 
a transient charring (or noncharring) conduction solution 
(providing component fluxes and surface equilibrium) . 

6. Edge boundary condition: 

(a) given total enthalpy, total pressure behind the shock, and 
pressure distribution around the body, 

(b) assigned boundary- layer-edge conditions including the possi- 
bility of an entropy layer. 

2. 3. 2. 2 Input And Output Data 

The input requirements for the BLIMP program are surprisingly few and 
simple considering the numerous options and the general applicability of the 
program. The options enumerated in Section 2. 3.2.1 are controlled by a single 
control card. Thermodynamic properties are provided by a thermochemical data 
deck (one card for each element and three cards for each molecular species) . 

A program is available for generating these data in the proper form. 
Multicomponent transport properties require only a set of diffusion factors, 

F i - Body shape is specified by nose radius and cone angle for sphere-cones 
or by r Q (s) for general axisymmetric bodies. The boundary-layer grid is 
established by a single set of rj-values plus the specification of a u/u e at 
one of the nodes (the a H (£) constraint) , which can be invariant from problem 
to problem as a consequence of the a TT parameter. In the absence of an en- 

ri 

tropy layer, boundary-layer edge conditions can be established from stagna- 
tion pressure, stagnation enthalpy, and pressure ratio around the body 
for the times of interest. For the case of surface mass and energy balances 
with surface equilibrium, no further input data is required with the excep- 
tion of material property data. In the other extreme of a completely specified 
(uncoupled) wall boundary condition, it is necessary to input T (or h ) 

and either 1) m and m or 2) and p v (or f ) , all for the 

y c '■v/ w w w 

times and streamwise stations of interest. 

The number of points across the boundary layer and the spacing between 
points, the number of streamwise stations and the distance between streamwise 
stations, and the number and type of elements (or base species) and candidate 
molecular and ionic species can all be selected arbitrarily. The size of 
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these arrays affect only dimension statements. It has been found (see Part 
III) that nominal 3-place accuracy can be obtained with about 7 spline points 
and 4-place accuracy with about 11 spline points if the spacing is optimized 
(e.g., T) of 0, 0.2, 0.5, 0.9, 1.4, 2.0, 2.7, 3.5, 4.4, 5.4, and6.5)and 
that poorly chosen nodal spacing can lead to convergence failures. 

The output from the BLIMP program consists of three parts: 1) a brief 

convergence history, listing a R , , damping factor, and maximum errors for 

each iteration; 2) data for current streamwise station (blowing rate, 
surface mass ablation rate, pyrolysis gas generation rate (or transpiration 
rate), mechanical removal rate, blowing parameters, elemental mass diffusive 
fluxes at the wall, wall heat flux, surface shear, skin friction coefficient, 
heat transfer coefficient, elemental mass transfer coefficients, momentum 
thickness, displacement thickness, shape factor, enthalpy thickness, elemental 
mass thicknesses, and Reynolds number per foot) ; and nodal data (distance 
from wall, r\, stream function, velocity ratio and its first and second 

derivatives with respect to T), total enthalpy and its first and second 

derivatives with respect to T], elemental mass fractions and their first and 

second derivatives with respect to r|,mole fractions of all molecular species, 
static enthalpy, temperature, density, viscosity, frozen specific heat, 
thermal conductivity, Prandtl number, molecular weight, and reference Schmidt 
number). A sample output is presented in Figure 3. 

2.3. 2.3 Storage Requirements and Computational Time 

The BLIMP program contains something in excess of 4000 instructions, 
including COMMON statements. The number of nodal points and the number of 
elements are the most critical dimensioned variables, in regard to both 
storage requirements and computational speed. To illustrate, the largest 
single block of numbers is the matrix of "nonlinear correction coefficients, " 
discussed previously, which is dimensioned [3 + (K + 2) n] by [3 + (K + 2 )n]. 
Furthermore, this matrix has to be inverted during each iteration. It is 
therefore significant that the solution procedure requires a relatively 
small number of spline points and usually converges in three or four itera- 
tions.* When dimensioned for seven nodal points, five elements plus elec- 
trons, 30 species, and 20 streamwise stations, the program fits on the 
Philco 212 without overlay but requires overlays on the IBM 7094 computer. 

By reducing the number of elements to four and molecular species to 25, it 
fits on the IBM 7094 without overlay but requires overlays on the smaller 
CDC 3200 computer 

The computational speed of the procedure is illustrated by the follow- 
ing examples. Equilibrium solutions for a sphere-cone configuration with 
coupled water injection into air involving 16 chemical species and seven 

*The relatively slow convergence in Fig. 3 ^(9 iterations) was a consequence 
of the 0.6 damping factor which was employed in this particular problem. 
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Figure 3. Typical Output from blimp Program. Stagnation Point Solution 
for Apollo Heat-Shield Material. Surface Equilibrium with 
Assigned Component Mass Fluxes 
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nodal points were obtained in approximately 2.5 minutes on an IBM 7094. These 
solutions included evaluation of edge conditions, a similar solution at the 
stagnation point, and nonsimilar solutions at ten additional stations. Problems 
with five elements and 30 species ' take nine seconds per iteration on the IBM 7094 
using overlays and 2 seconds per iteration on the Philco 212 without overlays. 

2.3.3 Sample Problem Solution 

The BLIMP program has been used to study a variety of problems. Several 
solutions are presented in Part III. A typical result is shown in Fig. 4. 
Profiles of velocity ratio, temperature, and shear function across a boundary 
layer into which a large quantity of Apollo heat-shield material is being in- 
jected are presented in Fig. 4(a). Mole fraction profiles are shown in Fig. 
4(b). These results were obtained for an assigned surface temperature and 
assigned component fluxes (m and m c ) and utilized a 30-component chemical 
model. A converged solution was obtained in seven iterations, starting with 
an air boundary-layer solution with the same wall temperature and same edge 
conditions but with no mass injection. 

2.4 SUMMARY 

The transformed nonsimilar laminar boundary-layer equations have been 
presented for a multicomponent boundary layer. These equations incorporate 
an approximate formulation for unequal diffusion and thermal diffusion coef- 
ficients for all species, a mixed equilibrium-nonequilibrium chemistry model 
and a one-dimensional model for radiation absorption and emission with angular 
dependent incident radiation at the boundary- layer edge. The formulation is 
suitable for coupling with a transient charring ablation solution and for 
matching with an entropy layer and a nonadiabatic inviscid flow field. 

A numerical procedure is then described for solving the above problem. 

Cubic spline functions are used to relate the velocity ratio, total enthalpy, 
and elemental mass fractions to the transverse boundary layer coordinate. 

The boundary-layer equations are integrated between neighboring nodes. This 
integration is primarily for algebraic convenience, the smoothness of the 
weighting function having been found to be relatively unimportant with regard 
to accuracy and convergence stability considerations. The streamwise deriva- 
tives (nonsimilar terms) are represented by conventional three-point finite- 
difference relations. 

A boundary- layer computer program, termed the boundary layer integral 
matrix procedure (BLIMP) , has been developed which utilizes this numerical 
procedure for solution of the nonsimilar, equilibrium, multicomponent, planar 
or axisymmetric, laminar boundary layer. This program applies to general 
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chemical systems and treats a variety of surface boundary conditions, indu- 
ing coupling to a transient charring ablation computational procedure. Accu- 
rate solutions have consistently been obtained with relatively few nodal 
points and convergence has been rapid, with the result that computational 
speed is a program virtue. 


SECTION 3 

ANALYSIS AND COMPUTATIONAL PROCEDURE FOR 
CHARRING MATERIAL RESPONSE 


3 . 1 INTRODUCTION 

Analysis of a complete transient ablation problem necessarily involves 
a computation of the internal thermal response of the ablating material. A 
substantial part of the total effort under the present program has been de- 
voted to the development of a computer code for in-depth response prediction 
of a charring material, suitable for coupling to the boundary layer program 
(BLIMP) described in Section 2 above. This section of the summary report 
describes the in-depth program (CMA) and related analysis. Section 5 below 
describes aspects of coupling procedures . 

3 . 2 PROBLEM DEFINITION 
3,2.1 General Description 

The basic problem is to predict the temperature and density histories of 
a thermally decomposing material exposed to some defined environment which 
supplies heat and which may chemically erode the material surface. 

The general prediction problem may conveniently be divided into two 
parts: the construction of a scheme for computing the in-depth behavior, and 

the specification of the heated surface boundary condition. The present report 
is mainly concerned with the first problem, although the second topic is also 
given extensive discussion. It may be noted in passing that for quasi-steady 
ablation problems (constant wall temperature, steady recession rate, invariant 
temperature profile with respect to the moving surface) , the details of the 
in-depth solution are not necessary for determining the surface temperature 
and the recession rate. The transient problem, on the other hand, does re- 
quire a complete in-depth solution, and hence is a much more elaborate prob- 
lem. 
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The physical problem may be illustrated as follows: 


heated surface 


char or residue 
pyrolysis zone 

virgin plastic 


As the material is heated, one or more components of the original composite 
virgin material pyrolyzes and yields a pyrolysis gas, which percolates away 
from the pyrolysis zone, and a porous residue, which for most materials of 
interest is a carbonaceous char, possibly reinforced with refractory fibers 
or cloth. 

Superimposed on this basic problem may be a number of even more complex 
events. The pyrolysis gases percolating through the char may undergo further 
chemical reactions among themselves and may react with the char, either 
eroding it or depositing additional residue upon it ("coking") . The char 
itself may collapse or fragment from mechanical or thermal stresses, and the 
refractory reinforcements may melt or suffer mechanical damage. Finally, 
various constituents of the residue structure may react chemically with each 
other, changing the nature of the char, and various mechanical forces may 
remove material from the surface. 

Despite these complexities, it is found that the "simple physics" des- 
cribed by 



virgin plastic — ► char + gas 

underlies a wide range of problems of technical interest, and for a great 
many materials, such as carbon phenolic, graphite phenolic, and wood, consti- 
tute all the events of interest. Such events as coking, mechanical erosion, 
melting, and subsurface reactions (other than pyrolysis) are less common 
and generally characterize specific problems. 

Therefore in any effort to compute the in-depth response of pyrolyzing 
materials the first order of business is to characterize the heat conduction 
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and the primary pyrolysis reaction which have useful generality. Particular 
details of special char chemical systems can then be superimposed upon this 
general computational scheme as required. The present effort has been mostly 
devoted to the general conduction pyrolysis problem. The numerical details 
will be described below. 


3.2.2 Differential Equations 

For the basic in-depth solution, it is assumed that thermal conduction 
is one-dimensional; however, the cross-section area (perpendicular to the 
conduction flux) is allowed to vary with depth in an arbitrary manner. This 
corresponds to a thermal stream tube. Furthermore, it is assumed that any 
pyrolysis gases formed are in thermal equilibrium with the char. Coking or 
further chemical erosion are not presently included in the computational pro- 
cedure and thus are excluded from the present discussion. An analysis has 
been developed for including these effects, however, and is discussed in Sec- 
tion 3.5 below. Thus, in the present discussion, it is assumed that the py- 
rolysis gases do not react chemically with the char in any way. Finally, any 
pyrolysis gas formed is assumed to pass immediately out through the char, 
that is, it has zero residence time in the char. Cracking or other chemical 
reactions involving only the pyrolysis gases may be simulated with an appro- 
priate gas enthalpy- temperature relation. 

The one-dimensional energy differential equation for this problem is 
readily formulated as 
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where p is the density, k is the thermal conductivity, A is area, h is 
enthalpy, and m^ the local gas flow rate. 

The conservation of mass equation is 




(36) 


Evaluation of this expression requires a specification of the decomposition 
rate 3p/d0)^. A great amount of laboratory pyrolysis data suggests that the 
decomposition rate may be taken as an Arrhenius type of expression 
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and for even greater generality it has been found useful and sufficient to 
consider up to three differently decomposing constituents 


p = np A + p b ) + (i - n p c 


( 38 ) 


where each component is governed by a relation of the form of Eq. (37) 



B . e 


-E i /RT 



i = A , B, C 


(39) 


For example (p, + p ) might be the density of resin (or analogous binder) 

1 A B 

in the ablating material, p c would be the density of the reinforcement and 
T the volume fraction of resin in the virgin plastic composite. 

It is possible to handle the decomposition in other ways than by Eq. (37) 
A popular simplification is to treat density as a function of temperature only 
An even more drastic simplification converts the virgin material to complete 
char at one particular "charring temperature." Other techniques specify some 
char thickness as a function of time or of heating rate. All of these sim- 
plifications are, of course, open to objection. Equation (37) is not only 
the most realistic physically, but is usually easy to handle in computation. 

3.2.3 Boundary Conditions 

Suitable boundary and initial conditions for the set of Eqs . (35) through 

(39) may be readily formulated. The boundary conditions at the front and back 
faces of the ablating material are usually surface energy balances. Of these, 
the front or "active" surface boundary condition is the most complex. It is 
handled in slightly different ways depending on which boundary layer treat- 
ment is being coupled to the in-depth response program. 

Basically, the surface energy balance may be pictured as 


i — 

w 

q rad j 

in 

, q rad , 

out 

, (pv) w h w 

Ev* 


■ ■ -J 

^cond 

i J 

m h 
c c 



m ti 

g g 


33 



where the indicated control volume is fixed to the receding surface (see Eq. 
(30)). Energy fluxes leaving the control volume include conduction into the 
material, radiation away from the surface, energy in any flow of condensed 
phase material such as liquid runoff, and gross blowing at the surface. En- 
ergy inputs to the control volume include radiation in from the boundary 
layer and enthalpy fluxes due to char and pyrolysis gas mass flow rates. The 
final input in the sketch is denoted — q^. It includes all diffusive energy 
fluxes from the gas -phase boundary layer (given by Eq. (18)) . If the in-depth 
response computation is being coupled to an exact boundary-layer solution, the 
term -q a is obtained from that solution procedure and is, of course, a 
complex function of the boundary-layer structure. If, on the other hand, 
the in-depth response is being coupled to a simplified boundary-layer scheme, 
such as a convective film coefficient model, then the term -q assumes the 
form of a correlation equation. Section 5 below contains a further discussion of 
this aspect of the total computation for these two approaches. 

For the present, it suffices to note that computation of the surface 
energy balance requires the following information from the in-depth solution: 

a. The instantaneous pyrolysis gas rate delivered from in-depth to 
the surface, it 

g 

b. A relation between the surface temperature and the rate of energy 

conduction into the material, a 

^cond 

With these two pieces of information the surface considerations allow deter- 
ation of char consumption rate m^ and surface temperature T . It 
will be useful to keep in mind that, from this point of view, the purpose of 
the in-depth solution at any instant is to provide information about m 

g 

and q cond^ T w^* In some circumstances, of course, it is of interest merely 
to specify the heated surface temperature and recession rate, in this case, 
no surface energy balance is required. 

It is usually of interest to have only one ablating surface. The back- 
wall or nonablating wall boundary condition may be modeled with a film coef- 
ficient heat transfer equation. 

3.3 FINITE DIFFERENCE SOLUTION PROCEDURE 

3.3.1 Introduction 

Section 3.2.2 above sets forth the governing differential equations 
whose solution is required to define the internal response of the charring 
material. As in many other problems, however, the differential equations 
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cannot be solved in general, and it is necessary instead to solve finite dif- 
ference equations* which model the differential equations and, the analyst 
hopes, retain the same mathematical properties as the original differential 
equations. A number of plausible difference equations can be proposed, and 
without the benefit of actual experience it is generally impossible to select 
any particular differencing scheme as superior to any other. In the past few 
years a few general differencing principles have been made reasonably clear, 
however, so that the analyst is not completely in the dark. The following 
section offers some background on this topic. 

3.3.2 Differencing Philosophy 

This section sets down the general principles upon which the finite dif- 
ferencing of the governing equations is based. These principles have proved 
sound and useful, particularly for complex problems. 

In common with all difference procedures, the area of interest (here, 
the charring material) is divided into a number of small zones, each consid- 
ered to be homogeneous. All derivatives in the governing differential equa- 
tions are then replaced by some difference expression from zone to zone. 

These zones, called nodes, thus provide the basic conceptual structure upon 
which the differencing procedure is based. 

The following principles of differencing and nodal sizing have been 
followed in the present programming effort; 

(1) The nodes have a fixed size. This avoids the slight additional com- 
putation complexity of shrinking nodes, and more importantly, makes principle 
(2) easier to satisfy, in addition to preserving a useful nodal spacing through- 
out the history of a given problem. 

(2) Since the nodes are fixed in size, not all of them can be retained 
if the surface of the material is receding due to chemical or mechanical 
erosion. From time to time a node must be dropped, and experience shows that 
it is much more preferable to drop nodes from the back (non-ablating) face 

of the material rather than from the front face. This means that the nodal 
network is "tied to the receding surface," and that material appears to be 
flowing through the nodes. This involves a transformation of differential 
equations (35) and (36) to a moving coordinate system and somewhat complicates 
the algebra of the difference equations modeled on these differential equations. 

*It is possible, of course, to use simpler schemes than finite difference 
equations. Integral analysis approaches, for example, have been tried. 

However, those techniques which have been employed have been of insufficient 
accuracy to be generally useful. 
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Disposing of nodes from the front surface, however, often leads to undesirable 
oscillations . 

(3) The difference forms of derivatives are kept simple and are formed 
so as to provide a direct physical analog of the differential event leading 
to the derivative. This approach may be contrasted to those approaches which 
seek elaborate difference approximations to derivative expressions. Experi- 
ence shows that the scheme advocated here, while sometimes at a minor dis- 
advantage in accuracy, greatly simplifies the attainment of a major objective: 
a difference scheme which conserves energy and mass. Many of the more elabo- 
rate difference schemes fail to meet these " simple" but crucial conservation 
criteria, and hence frequently converge to erroneous or spurious solutions. 

(4) The difference equation for energy is formulated in such a way 
that it reduces to the difference equation for mass conservation when tem- 
peratures and enthalpies are uniform. Any lack of consistency between the 
energy and mass equations complicates, and may entirely defeat, convergence 
to a meaningful result. 

(5) The difference energy equations are written to be " implicit" in 
temperature. That is, all temperatures appearing are taken to be "new" 
unknown temperatures applicable at the end of the current time step. It is 
well established that implicit procedures are generally more economical than 
explicit procedures, at least for the majority of ablation problems of inter- 
est in the current work. 

(6) In contrast to point (5) , the decomposition relations are written 
as "explicit" in temperature. To implicitize temperature in these highly 
nonlinear equations necessarily involves either a time-consuming iteration 
procedure or an elaborate linearization. 

(7) Since experience has shown that material decomposition rates are 
strongly dependent on temperature, it is highly desirable to perform the 
mass balance operations in a different, tighter network than that used for 
the energy balance equations. For greatest generality and utility, the num- 
ber of these mass balance "nodelets" per energy balance "node" are 

freely selectable. 

3.3.3 Array of Difference Equations 

The actual derivation of the necessary finite difference equations is a 
complicated and tedious mass of detail best left to the detailed development 
given in Part II of the present series. For this summary report, only the 
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final result of the energy equation differencing is of interest, since the cou- 
pling procedure description given in Section 5 below requires an understand- 
ing of this step. 

With an implicit temperature formulation, it may readily be seen from the 
energy differential equation (35) that the difference form of the energy equa- 
tion for a given node at any time step will involve the "new" unknown tempera- 
ture of that node and the two adjacent nodes. The equation for the last node, 
however, will have only two unkowns, since the adjacent temperature for that 
node is the known reservoir temperature. Similarly, the first node equation 
involves only two unknown temperatures, but it also includes the unknown quan- 
tity q which ultimately will be determined by the surface energy balance. 

If we arrange all the energy relations in order we obtain an array of 
the form 


B T* + C T' 
11 i < 


^cond 


AT' + B T' + C T* 

ax 2 a 2 3 


D 
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AT’ + B T’ + C T* 

3 2 3 3 3 4 
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AT 1 + 

4 3 


B T' + C T 1 

4 4 4 5 


D 


4 


Vi-, + Vn + ° = f(T res> 

(40) 

where the coefficients B N> C^, and are given by complicated 

algebraic expressions involving nodal thickness, thermal conductivities, gas 
flow rates, and old temperatures. 

3.3.4 Reduction of Array 

It is now possible to see clearly what needs to be done for each time 
step Ae of the solution in order to prepare for coupling to the surface 
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energy balance. First, using the current values of p n , surface recession rate, 

and T n , the mass relations (38) and (39) can be solved for each node n, 

yielding "new" gas flow rates m . 

^n 

This information may then be used to compute the coefficients of the 
tri-diagonal energy equation matrix. Once this matrix is set up, the required 
surface energy relation q cond s ^cond^w^ ma ^ be obtained directly, as des- 
cribed in the next section. 

The next step in the solution procedure is to eliminate one unknown tem- 
perature from each equation in the array (Eq. (40)). This can be done by 
eliminating T^ from the next to last equation and proceeding sequentially 
upward to the top equation, eliminating the highest- indexed unknown tempera- 
ture from each equation of Set (40) in turn. The * resulting reduced set has 
the form 


B*T ' 

l i 


^cond 


AT’ + B*T' = n* 

3 12 3 3 


AT* + B*T' = D* 

3 3 3 3 3 


AT’ + B*T ' 

4 3 4 4 


D* 

4 


Vb-i + B N T N = F 4 < T res> ^ 41 ) 

3.3.5 Coupling to Surface Energy Balance and Final Step 

The top equation of this set relates q to T 1 , where T’ is the 

cone* ^ ^ 

new surface temperature. This is the required relation for use in the sur- 
face energy balance coupling to whatever boundary-layer model is being used, 
as discussed in Section 5 below. This coupling operation balances the 
surface energy terms and thus determines tie new surface temperature T* . 

Since T x * is now known, the second equation of Set (41) yields T^ directly, 
then the third equation yields , and so on until the new temperature set 
is complete. 
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As a final step, new values for temperature dependent properties can be 
selected for each node and the entire system is then ready for a new time 
step, beginning with the decomposition event. 

3.4 CHARRING MATERIAL ABLATION (CMA) PROGRAM 

The Charring Material Ablation Program is a coded procedure for calcu- 
lating the in-depth thermal response of a charring, ablating material. The 
solution is obtained through the difference equations discussed in the pre- 
ceding subsections. The program is described briefly in Section 3.4.1 and 
some sample problem solutions are presented in Section 3.4.2. Complete de- 
ceptions, user's manual, and flow charts are given in Refs. 12, 13, and 14. 

3.4.1 Program Description 

3. 4. 1.1 Program Objectives 

The program produces in-depth temperature and density histories, plus 
surface recession rate as a function of time. In addition to this basic out- 
put, the program outputs a number of integrated energy terms and various mate- 
rial property data of interest. Section 3. 4. 1.5 below gives a more detailed 
description of the program output. 

3. 4. 1.2 Program Capabilities 

The Charring Material Ablation Program is an implicit, finite-difference 
computational procedure for computing the one-dimensional transient transport 
of thermal energy in a three-dimensional isotropic material which can ablate 
from a front surface and which can decompose in depth. Decomposition reac- 
tions are based on a three-component model. The program permits up to eight 
different backup materials of arbitrary thickness. The back wall of the com- 
posite material may transfer energy by convection and radiation. 

In one program configuration, the ablating surface boundary condition 
may take one of three forms: 

Option 1 — Film coefficient model convection— radiation heating with 
coupled mass transfer, including the effects of unequal 
heat and mass transfer coefficients (non-unity Lewis 
number) and unequal mass diffusion coefficients. Surface 
thermochemistry computations need not presume chemical 
equilibrium at the surface, and may allow for melting and 
liquid phase removal at the surface. Chemical state programs 
for providing this thermochemical surface boundary condition 
are discussed in Section 4.4. 
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Option 2 - Specified surface temperature and surface recession rate 

Option 3 - Specified radiation view factor and incident radiation 

flux, as functions of time, for a stationary surface. 

Any combination of the first three options may be used for a single 
computation. Option 3 is appropriate to cooldown after termination of con- 
vective heat input and is often useful in conjunction with Options 1 and 2. 

In another configuration, the program may be coupled to the boundary 
layer integral matrix procedure (BLIMP) program. In this arrangement, the 
total assembly is designated the CABLE program and is described in Section 5 
below. 

The program permits the specification of a number of geometries. in the 
most general case, area may vary arbitrarily with depth. Special cases in- 
clude: 

(1) Plane 

(2) Cylindrical or annular, with heated surface either inner or outer 

(3) Spherical or spherical shell, with heated surface either inner 
or outer. 

The rear surface of the last node may be specified as insulated, or may 
experience convective and radiative heat transfer to a "reservoir" at a spe- 
cified reservoir temperature if a rear surface convection coefficient and an 
emissivity are input. 

Material properties such as thermal conductivity, specific heat, and 
emissivity are input as functions of temperature for virgin plastic and char. 
For partially decomposed material, the program performs an appropriate averag- 
ing to determine effective material properties. 

3.4.1. 3 Solution Procedure 

The basic solution procedure is by the finite difference approach dis- 
cussed in Section 3.3. For each time step, the decomposition relations are 
solved and then the in-depth energy fluxes constructed in general terms. 

These are then harmonized with a surface energy balance (if a surface energy 
balance option is being used) and the in-depth temperatures determined. New 
material property values are set up and the solution is ready for the next 
time increment. 

3. 4. 1.4 Output Information 

The CMA program outputs instantaneous mass ablation rates and blowing 
parameters for char and pyrolysis gas, total integrated mass ablation of char 
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and pyrolysis gas, total recession and recession rates of surface, of the char 
line, and of the pyrolysis line. It also outputs the surface energy flux 
terms, namely, the energy convected in, energy radiated in, energy reradiated 
out, chemical generation, and conduction away (c[ cond ) - Further, it describes 
how the input energy of q cond i s "accommodated" or "partitioned" in the 
solid material. Part of the energy is consumed in decomposing the plastic, 
part is consumed in sensible enthalpy changes of the solid, and part is 
"picked-up" by the pyrolysis gases as they pass through the car. Thermocouple 
and isotherm output can also be called for. A typical output is presented in 
Figure 5. 

3 .4 . 1 . 5 Storage Requirements and Computational Time 

The storage requirements for the CMA program depend strongly upon the 
coupling mode in use. Coupling to a film coefficient model for the surface 
energy balance (option 1) involves considerable table storage such that the pro- 
gram will barely fit a 32,000-word machine with full table sizes. In certain 
cases a reduction in table sizes will allow the program to fit on a smaller 
machine. As a subroutine to the CABLE program or use of Option 2 or Option 3 
eliminates the need for storing extensive boundary condition tables. In these 
cases, the CMA program requires less than 8000 words of storage. 

Computation time depends, of course, on the problem being computed, but 
experience to date indicates that CMA computations run in roughly 1/3 of real 
time for "typical" charring material problems, for machines of the IBM 7094 
speed class. 

3.4.2 Sample Problem Solutions 

As an illustration of the general performance of the charring material 
computer program. Figure 6 (a) presents a graphic representation of the in-depth 
density history for a nylon— phenolic material exposed to a typical reentry 
environment. Figure 6(b) shows some in-depth thermocouple temperature response 
predictions for the same problem. Figure 7 , from a different problem, shows 
a machine made plot generated by a plot routine coupled to the CMA program. 

3.4.3 Concluding Remarks for the Charring Material Ablation Program 

The preceding sections have described the analysis and an associated 
computer program for the calculation of the in-depth response of a charring 
or pyrolyzing material. The general objective of the development effort has 
been to produce a computation scheme which accounts for those physical events 
common to a wide range of technically important applications, so that the re- 
sulting program has as much generality and flexibility as possible. To this 
end, the analysis accounts for the basic in-depth pyrolysis problem, which 
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Figure 5. Typical Output from CMA Program. Stagnation Point Solution for 
Apollo Heat-Shield Material During SA 202 Trajectory 
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is observed to be common to a wide range of problems, and excludes coking 
(char densification) , thermal expansion, condensed phase char rate-controlled 
reactions, and mechanical damage mechanisms. All of these are specific to 
particular materials or material types. For such materials, the basic pro- 
gram can be modified to include these special effects. 

The basic program generates a one-dimensional in-depth solution, but the 
cross-sectional area of the material analyzed may vary with depth (thermal 
stream tube) . Pyrolysis may occur through three distinct Arrhenius-type 
kinetic reactions. 

An important feature of the program is the range of physically realistic 
boundary conditions available for the heated surface. These include 

(1) Specified temperature and recession rate 

(2) Radiation energy balance with zero recession and no convection 
(cool down or soak out) 

(3) Coupling through a film coefficient model to surface thermo- 
chemistry solution including general heterogeneous equilibrium, 
or heterogeneous equilibrium modified by certain rate controlled 
reactions, both models including the effects of the melting and 
total removal of surface species formed at temperatures above 
their melt or fail temperatures 

(4) Coupling to a general, nonsimilar, multicomponent boundary-layer 
solution including heterogeneous kinetic effects, with surface 
melting or failing. 

This range of possibilities offers opportunities for economy during routine 
in-depth studies or during computations for which film coefficient models 
are adequate, while preserving the capability of doing very accurate coupled, 
simultaneous boundary layer and in-depth solutions. 

Other features of particular importance which are worth stressing may be 
enumerated as follows: 

Simple difference formulation . The difference equations employed to 
effect the in-depth solution are kept simple and are designed to conserve 
energy and mass without fail. This helps avoid convergence to spurious 
solutions . 

Implicit temperature formulation for energy, explicit temperature for - 
mulation for decomposition . This approach keeps the energy solution stable 
and economical and avoids the complexities of a totally implicit solution. 
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Three- component Arrhenius expression for pyrolysis kinetics . This 
expression is believed sufficiently general to accurately represent the 
degradation kinetics of most ablation materials of interest. 

Freely selectable pyrolysis gas kinetics model . The history of the 
pyrolysis gas as it passes through the char is a key aspect of any solution. 

The analysis scheme used here does not involve any particular assumption 
about the pyrolysis gas history, but employs only an input temperature-enthalpy 
relationship. This relationship may be derived from any model the user feels 
applicable . 

Final conclusions about the general accuracy and utility of the CMA in-depth 
routine will be deferred until Section 5.2.4 below where the program can be con- 
sidered as coupled to a thermochemical boundary condition. For the present it may 
merely be noted that the in-depth aspects alone of the solution have been sub- 
jected to extensive scrutiny for cases of known surface temperature and recession 
rate. Agreement with both analytical solutions and with in-depth thermocouple 
data has generally been excellent, although, of course, when dealing with charring 
ablators, adjustments in thermal conductivity and pyrolysis gas enthalpy data are 
frequently necessary to force agreement between predictions and thermocouple data. 
These two parameters have such powerful effects on the solution that the thermo- 
couple data must be supplemented by in-depth observations to detect any omissions 
in the analytical model. The chief possibilities here are coking and reinforcement- 
char chemical reactions; both effects have occasionally been observed. The con- 
clusion remains, however, that the program accurately models the fundamental 
pyrolysis and energy transport events it was intended to model. 

In conclusion, the in-depth analysis presented here and programmed as the CMA 
program has been applied to a wide range of materials of technical interest with 
excellent results. The program appears to be thoroughly checked out and fully 
operational . 


3.5 SUBSURFACE COKING REACTIONS 

The CMA program described above is a mathematical analog of the surface 
and subsurface thermochemical events which are common to most reinforced 
organic ablation materials. These events consist basically of (1) energy 
and mass transfer in a material which experiences subsurface decomposition 
of an organic polymer into a pyrolysis gas and a char residue, and (2) the 
energy, mass, and chemical species transfer at the ablating surface which 
dictate the magnitude of surface recession. As discussed earlier, additional 
complications may be of importance when consideration is given to certain 
particular types of ablation materials. These complications may include such 
events as subsurface reactions between reinforcing fibers and the carbonaceous 
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char, mechanical failure of the char layer resulting from pressure gradient 
and thermally induced stresses, and the thermal cracking of pyrolysis gas 
products resulting in carbon deposition in the char layer beneath the heated 
surface. Phenomenological models have been postulated and some experimental 
data has been obtained for both subsurface fiber— carbon interactions (Ref. 15) 
and mechanical failure of the char layer (Ref. 16) • Little attention has 
been directed toward the construction of a phenomenological model to repre- 
sent subsurface pyrolysis gas cracking with attendant carbon precipitation 
upon the char layer (coking) . 

An analysis has been conducted and finite difference relations suitable 
for computer coding have been developed for representing the response 
of a charring ablation material which may undergo subsurface "coking" of the 
pyrolysis gas. Coking reactions, as employed here, refer to the precipitation 
of carbon from the hydrocarbon— containing gaseous pyrolysis products with 
attendant deposition upon the char, and the reverse reaction evidenced by 
erosion of the carbonaceous char accompanied by the addition of carbon to the 
gaseous pyrolysis products. Forward and reverse coking reactions are of 
interest because the permeability of the char layer is decreased by coking 
which may result in high gas pressure in depth. High pressure in depth may 
give rise to excessive char stress which may produce catastrophic failure 
of the char layer. The technique considers internal pressure build-up result- 
ing from pyrolysis product flow through a char layer of variable permeability. 

Details of the development and finite difference formulation are presented 
in Part VI of the present series. A brief description of the adopted phenomeno- 
logical model and some comments regarding the finite difference formulation are 
given in the following two subsections. 

3.5.1 Description of the Physical Process 

The physical process to be characterized represents an extension to that 
represented in the CMA program described above, in that after the pyrolysis 
gas is formed, further mass transfer between the pyrolysis gas and char may 
occur. The resultant change in pyrolysis gas and char composition is evaluated 
and the effect of this change upon subsequent events, both below and at the 
heated surface are considered. The generalized model selected to represent the 
charring material is described first. This is followed by a description of the 
type of reactions that may be considered and the rate equations selected to 
represent the rate at which these reactions may proceed. 

In its undecomposed state the ablation material is taken to be composed 
of two basic types of constituents: (1) inert, and (2) reactive. The inert 

constituents will consist of materials which are not permitted to undergo 
molecular changes in depth, e.g., silica or other metal oxide reinforcements. 
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The reactive constituents may consist of organic materials, carbon or graphite 
reinforcements, and water of crystalization of reinforcing fibers, for example. 
Carbon and graphite are included in the list of reactive constituents because 
they may be vaporized in depth or be eroded chemically by the gaseous products 
of other reactive constituent pyrolysis products. The following, idealized, 
irreversible reaction characterizes the initial decomposition of the composite. 

Inert + Reactive -► Inert + Carbon (S) + Gas (0) (42) 

As noted from this relation, the inert constituent does not take part in the 
reaction, but is simply transported from a constituent in the virgin plastic 
to a constituent in the initial decomposed material. The reactive consti- 
tuents, on the other hand, do undergo a change in molecular configuration and 
phase; however, the products of this initial reaction may consist of only two 
constituents, solid carbon, and an initial pyrolysis gas (gas(O)). The reac- 
tion should be looked upon as one which splits the virgin material into three 
distinct parts, each having a fixed quantity of chemical elements. 

I'he initial off-gas elemental composition may then be obtained by sub- 
tracting the quantity of chemical elements contained in a laboratory-produced 
char from the chemical elements contained in the virgin plastic. After the 
initial decomposition gas is formed, it will percolate through the char layer 
toward the heated surface . This will result in an increased gas temperature 
and decreased pressure. The change in pressure and temperature will cause the 
initial gas products (gas (o) ) to undergo numerous chemical reactions as they 
pass through the char. The reactions considered fall into three general cate- 
gories ; 

1. Decomposition of the gas including thermal decomposition of high- 
molecular-weight hydrocarbons and dissociation of C0 9 , H 0 0, and H 0 
for example. 

2. Further decomposition of the hydrocarbons resulting in precipitation 
of carbon (coking) on the adjacent char passages resulting in char 
density buildup. 

3. Chemical erosion of the char layer (below the heated surface) by 
the gases including carbon vaporization resulting in a char density 
reduction near the heated surface. 

The three reaction regimes are represented in the following sketch for decom- 
position of a hypothetical material. 
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Char 

Erosion 



The classification of reactions above corresponds to the order in which the 
various types of reactions would be expected to occur as the gas passes through 
the char. The first class of reactions may be looked upon simply as the gas 
reacting with itself so no change in the concentration of chemical elements 
in the gas results. The second two classes of reactions, however, will result 
in a transfer of carbon elements between the gas and the porous char. With 
regard to the elemental composition change these reactions may be considered 
reversible . 

2 

gas (0) gas + carbon (S) (43) 

3 

where the forward reaction corresponds to moderate temperature, type 2 reac- 
tions, and results in a precipitation of carbon from the initial pyrolysis 
gas. The subsequent, type 3, high temperature reactions result in char ero- 
sion with attendant addition of carbon to the gas. 

An understanding of the detailed kinetic mechanisms required to charac- 
terize these reactions is not presently in hand; however, some qualitative 
information is available upon which a crude model may be formulated. The 
specific information relating to^each reaction type is presented briefly and 
the physical model adapted to characterize each reaction regime is described 
in the following paragraphs. 

Type 1 Reactions - If the subsurface composition is computed on the basis 
of chemical equilibrium considerations, a far more dense char is predicted to 
occur than is observed from char density measurements (Ref. 17) • It may be 
concluded either that condensed phase carbon is formed but does not stick to 
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the char, or that the high molecular weight hydrocarbons do not decompose 
according to the dictates of chemical equilibrium. The latter possibility 
seems more probable but conclusive experimental evidence on this matter is 
lacking. 

Type 2 Reactions - A certain amount of experimental evidence exists 
(Refs. 18 and 19) which indicates that char densification may occur between 
the organic decomposition zone and the heated surface. It appears that this 
char densification is a result of deposition of condensed phase carbon from 
the hydrocarbons in the organic pyrolysis gas products. This seems reasonable 
since, as indicated above, the gases contain far more carbon than would exist 
if equilibrium were achieved. As the gases approach the heated surface their 
temperature is increased and the rate at which equilibrium is approached 
increases. In the present study, the coking rate is expressed as the product 
of a forward rate coefficient and a carbon mass fraction "coking potential." 


coke 


V K cg - K cgE> 


(44) 


where the forward rate coefficient is expressed in Arrhenius form, and the 
driving potential is represented by the difference between the elemental 
carbon mass fraction of the gas and that which would exist if equilibrium 
were achieved. The coking rate equation has the following features: 

1. It is simple enough to be included practically in a charring 
ablation solution. 

2. It approaches the coking rates which would be predicted by more 
detailed complete equations when kinetics are relatively slow or 
fast. 

3. It is based upon parameters vdiich may be controlled in a laboratory 
experiment to derive data on the coking process. 

As noted from the postulated "coking" rate equation, the coking rate will 
approach zero as the pyrolysis gas approaches chemical equilibrium with the 
char. This marks the onset of regime 3 which is characterized by addition 
of carbon to the gas from the char. 

Type 3 Reactions - in the event local char layer temperatures much in 
excess of 4500° R are achieved, it is probable that chemical equilibrium will 
be achieved between the pyrolysis gas and the char, in which case the coking 
potential will reach zero. For this reason type 3 reactions are presumed to 
occur in chemical equilibrium and sufficient carbon will be added to the gas 
from the char to maintain this equilibrium. This char erosion will result in 
a char density decrease which may, in an extreme case, cause the char density 
to approach zero. 
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3,5.2 Solution Procedure 

The differential and finite difference formulation of the energy, mass, 
and momentum transfer events associated with the above phenomenological model 
is accomplished in much the same manner as with the CMA program described 
earlier. The primary feature which distinguishes between the two formulations 
is associated with the treatment of the subsurface species conservation equa- 
tions and evaluation of the char material state. The species conservation 
equations include evaluation of the pyrolysis gas elemental carbon content 
as results from the integrated effect of coking events. 

In order that the energy exchanges associated with coking reactions be 
properly evaluated it is necessary to accurately assess the thermodynamic 
state of both the pyrolysis gas and the char material. This is accomplished 
by expressing the pyrolysis gas enthalpy as a function of pressure, tempera- 
ture, and elemental carbon content. The char state is expressed in terms of 
temperature and the relative quantities of each constituent; inert, reactive, 
and carbon. Because the ablation material composition depends upon the rela- 
tive quantity of each of three constituents it is not possible to express the 
material state in terms of temperature and density alone, as in the CMA pro- 
gram. Because of this, it is most convenient to express both the energy and 
mass conservation equations in terms of nodal mass rather than nodal density. 

In addition to the inclusion of energy and mass transfer events associated 
with subsurface coking relations, an empirical expression is included for 
evaluating the pressure distribution through the char layer as results from 
pyrolysis gas flow through the char layer of variable permeability. 

3.6 CONCLUDING REMARKS 

The above sections have described the basic physical problems of the 
charring ablator and have given the governing differential equations . The 
finite difference solution procedure for solving for the in-depth response 
of a charring ablator was sketched out. The manner in which this in-depth 
solution could be coupled to either a film coefficient boundary-layer model 
or to a complete boundary-layer solution was indicated; this subject will be 
dealt with in greater detail in Section 5 below. Sample solutions were cited 
illustrating the general applicability of the analysis in its computer program 
form. 

A supplementary analysis was described which broadened the physical prob- 
lem to include the effects of carbon deposition from the pyrolysis gases 
(coking) and the reverse char erosion effect. The resulting solution proce- 
dure differs in detail from the no-coking procedure, since a new physical quan- 
tity (carbon content) must be traced, but the general level of problem complexity 
is not appreciably higher than in the no— coking problem. The coking analysis 
has not been programmed for machine computation. 
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SECTION 4 


ANALYSIS AND COMPUTATIONAL PROCEDURES 
FOR EVALUATION OF CHEMICAL STATE 


4.1 INTRODUCTION 

In the study of high energy boundary layer phenomena, thermochemical pro- 
cesses can be of dominant importance. This is particularly true when these 
boundary layers interact with chemically active surfaces, in the present 
study, interest is directed toward the prediction of thermochemical response 
of a heat shield during superorbital reentry. The requirements for evaluating 
the chemical state of homogeneous and heterogeneous systems in this problem 
are extensive. These requirements include the determination of the chemical 
state after normal or oblique shock wave compression, during the isentropic 
expansion of the inviscid shock layer gases, within the boundary layer, and at 
the chemically active surface. In the last two instances, these state calcu- 
lations are coupled with complex mass balance relations. Many chemical state 
solution procedures have been documented to treat reasonably standard closed 
systems, such as those associated with expansion processes. For open systems 
only a few direct solution procedures have been documented. Because of the 
number of requirements imposed upon the chemical state routines in the present 
study, the general treatment of a variety of chemical systems became a major 
effort. The inclusion of a general kinetic model, ionization, and the exten- 
sive bookkeeping associated with the downstream introduction of new species, 
is of major importance in the formulation of the general problem necessary 
for thoroughly treating the coupled boundary layer problem. 

The chemical state procedures adopted as a part of this effort are de- 
scribed in Part V of this series of reports and summarized in the following 
sections. The basic relations are presented in Section 4.2, whereas in Sec- 
tion 4.3 the solution procedure is discussed. These techniques have been 
built in greater or lesser extent into the equilibrium surface thermochemistry 
(EST) program, the Aerotherm chemical equilibrium (ACE) program and certain 
special modifications of it, and the chemical state subroutines to the BLIMP 
program. Section 4.4 specifies the status of these routines and the extent 
to which the general formulation presented herein has been implemented. In 
brief, the procedures are presently limited to equilibrium, except that se- 
lected species can be considered as frozen across the boundary layer and to 
undergo rate-controlled surface-catalyzed reactions or reactions with the 
surface material. 

4 , 2 ANALYTICAL APPROACH 

In this section, the analytical approach utilized to specify the chemical 
state of a system are summarized. Basically four types of relations can be 
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considered in a general open system. These are the equilibrium relations 
applying to those reactions which can be considered as generally equilibrated, 
the nonequilibrium relations for those reactions which can be (but are not 
necessarily) out of equilibrium, the mass balance relations, and those addi- 
tional state constraints imposed on the system. 


4.2.1 Equilibrium Relations - Totally Equilibrated Systems 

In a chemical system there will exist, in the general case, a set of in- 
dependent equilibrium reactions. All other equilibrium reactions will be 
equivalent, both physically and mathematically, to this independent set. It 
can be shown that in a completely equilibrated system the number of independent 
equations is usually equal to the number of molecules less the number of ele- 
ments. The modification of these relations for systems that are not completely 
equilibrated will be considered in Section 4.2.2. 

The selection of this set of independent reactions can be done arbitrarily, 
but it is convenient to establish some consistent technique. Most such tech- 
niques are based on the pre-selection of a set of species usually equal in 
number to the number of elements. The formation reactions of all other species 
from this base set represent the independent set of equilibrium reactions. The 
base species must be selected in such a fashion that no reaction can be written 
wherein reactants and products are all base species. Thus in the 0,H system, 

HO and H 2 02 represent an invalid base set whereas HO and 0, HO and H, etc., 
represent valid sets. It has been reasonably common practice to select the 
monatomic gases as base species, since the formulation of the formation reac- 
tions is particularly convenient. There are advantages, however, in selecting 
a more general set, particularly when chemical kinetics are important. Con- 
sidering a set of base species N^, formation reactions for the remaining Nj 
species are of the form 


Z v ji N i - N j 


where the 
Mathematically, 
of element k 


are the stoichiometric coefficients of the formation reactions. 

the v . . are obtained implicitly from the c v . (the atoms 
j 1 

in molecule j) by 


Z c ki v j± = c kj 


The set of independent formation reactions (Eq. (45)) for j ranging 
from ^ to N g , where is the number of base species and N g is the 

total number of species, can be used to formulate a set of equilibrium 
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constraints. At equilibrium, the second law requires that these independent 
reactions occur without change in free energy. Therefore 


= Y"* v . . G. 

Z_i 311 


(47) 


where the Gj are the partial molar free energies of the species, it is 
shown in Part V of this series of reports that equilibrium constant relations 
can be expressed as 


in K 



in Pj 


> v . . in p . 
3 1 *1 


(48) 


where p^ is the partial pressure of the j th species, T is temperature, 
R is the universal gas constant, the standard state free energy change of 
the formation reaction for species j is defined by 


AG° 


G 


o 

j 



(49) 


and the partial pressure of condensed species is taken as one atmosphere. 
The standard state free energy is a function of temperature only and is ob- 
tained for each molecular species from 



(50) 


where enthalpies are obtained relative to some chemical base state, often the 
elements in their most natural form at 298° K and one atmosphere (JANAF base 
state) . 

The stationary condition of the free energy at equilibrium expressed in 
Eq. (47) is consistent with the minimum free energy statement often utilized 
in seeking the equilibrium state. Although the formulation followed here 
differs from those followed in free energy minimization approaches, the ulti- 
mate numerics can reduce to an identical iterative solution procedure. 

The solution of the set of algebraic equations (Eq. (48) ) must be con- 
sidered in conjunction with other constraints including the pressure balance 

Pj = P (51) 

j 
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where the summation is over all gas phase species* The detailed solution pro- 
cedure will be considered only after all required relations have been discussed. 


4.2.2 Mixed Equilibrium-Nonequilibrium Relations 

When some reactions fail to equilibrate it is necessary to approach the 
selection of the independent set of equilibrium reactions with greater caution. 
In a general chemical system, certain sets of molecules can be treated as 
always equilibrated. Between these sets certain independent equilibrium and 
kinetically controlled interchange reactions may exist, A procedure for treat- 
ing mixed equilibrium- nonequilibrium systems is presented in this section. 

The following rules are established in order to organize the logic: 

1. Every species is assigned to one and only one set, 

2. A set may contain as few as one species. 

3. Each set has its own base species, i.e., that minimum number of 
species from which all other members of the set may be formed. 

4. Within each set all possible reactions between member species 
are equilibrated. 

5. Equilibrium interchange reactions involve species from more 
than one set. 


Consider, for example, eight species of the 0-H system: 0, H^O; H; H 2 ; 
°2 °3 ; H 2°2 w ^ ere fi- ve s ets are divided by semicolons. For these sets, 

the following base species are appropriate: 0, H 2 0; H: H 2 ; 0 2 ; HO where only 
the first set requires more than one base species. At this juncture only 
two independent equilibrium reactions have been formulated, namely 


1.5 0 2 


2 HO ^ H 2 0 2 


} (52) 


Two independent equilibrium interchange reactions might be included in this 
system, for example 

H + OH 4± H 2 + 0 

\ (53) 

0 + OH ** 0 2 + H ] 


The effect of these reactions is to reduce the number of base species by two. 
For example H 2 and 0 2 can be deleted. The remaining base species and the 
array of formation reactions coefficients, v jj_ # are therefore 
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1 

2 
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in the general relation of Eq. (45) . 

If it is assumed that all other interchange reactions are frozen the 
formulation of the equilibrium-nonequilibrium aspects of the program are con- 
plete. In the totally equilibrated chemical system, conservational constraints 
are often applied to the elements. in the system just presented, however, 
additional conservational constraints are required. In general these con- 
straints take the form 



3 


v . . n . 
3 


= a . 


(54) 


where n^ is the number of moles of species j in a unit mass of system and 
is a conserved variable relating to the '•elemental* 1 composition of a unit 
mass of the system. One such constraint is applied for each base species. 

In effect, the base species become the "elements" of the system, and their 
total masses can be treated as the conserved variable of the system. This 
generalized concept of the conserved "elements"* of the system is extremely 
important to the present development. 


In the general nonequilibrium system, certain kinetically controlled 
reactions will be important. For example, in the H-0 system 


H + H + M H 2 + M 




0 + 0 + M 0 2 + M 


} (55) 


H 4- H + 0 -► H 2 0 




The term "element" (in quotes) is used to refer to those atoms or groupings 
of atoms (i.e., grouped according to the base species formulae) which accord- 
ing to the equilibrium relations are conserved. 


56 




are three reactions of possible interest, M being any third body. The rates 
of these reactions can be related to the partial pressures of the reactants 
and products, the ultimate equilibrium constraint appropriate to the reaction, 
and the kinetic coefficient. With the general kinetic reaction in the form 


Z u ?» N j - Z 4»“i 

j j 


(56) 


its rate, R^, 


can be expressed generally by 



where k is the forward rate constant. The net effect of these reactions 
F ra 

is the modification of the "elemental'’ makeup of the system. The kinetic 
reactions cause a net increase rate (moles per unit volume) 


r. = ) ) (m P -h R )v..R (58) 

i L. J 3 m 3 m 31 m 

m j 

of "element" i. It is this relation which is introduced into the conserva- 
tional equations in order to establish the local state of the reacting chemi- 
cal system. 

The forward rate constant k™ , hopefully based on experimental data, 

m 

is represented with an Arrhenius type function 

k = B m exp (E a /RT) (59) 

m m 

v/here the exponential factor establishes the probability of a collision having 
energy in excess of the activation energy E a ^ and the factor B m represents 
a multitude of phenomena associated with the probability of success of a single 
collision (e.g., collision orientation). 

When kinetically controlled reactions approach equilibrium, difficulty 
is often encountered in the treatment of the relevant conservational equations. 
To understand the nature of this difficulty, and thus the means of avoiding 
it, it is instructive to consider the simple time dependent character of the 
H-0 system previously described. Recalling that represents the moles of 
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"element" i in a unit mass, and that r i represents the rate of production 
of moles of "element" i per unit volume, it follows that 

dou r^RT 

~Se~ = ™ p vT (60) 


At this point it is necessary to introduce another new concept. From 
the base species, a subset of base-base species can be obtained much 

as if all specified kinetic reactions were permitted to equilibrate. For 
this example 0 and H will be selected for this honor and the "formation 
reaction" for the remaining base species written as 


2H + 0 — ► H 2 0 
H + 0 — ► OH 


(61) 


or, more generally. 


I 




N. 


1 


(62) 


The reactions (Eqs. (61)) equilibrate if the third of the kinetic reactions 
of Eq. (55) and either of the other two reactions have infinite rates. As 
shown in Part V of this series of reports, of the Eqs. (60) can be re- 

placed by 



(63) 


where 



(64) 


For each base species i which is also a base-base species k, an 
equation of the form of Eq. (63) replaces the corresponding one of the form 
of Eq. (60)_. The other Eqs. (60) are maintained unaltered in the system of 
equations and still contain the kinetic expressions. These are referred to 
as the reactive mass balance equations. In the general case, a given reaction, 
m, will affect more than one of these equations. The consequences of this 
when such a reaction approaches equilibrium are discussed in Part V of this 
series of reports. By combining equations in a manner based on the a priori 
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selection of a set of controlling reactions, these consequences can be cir- 
cumvented in an effective fashion, 

4,2.3 Mass Balance Relations 

In the preceding sections the equilibrium and nonequilibrium relations 
have been summarized for a general chemical state. These relations are in 
themselves insufficient until other relations, in particular the mass balance 
relations, are imposed. In the case of kinetic control, the time dependence 
of the system must be equated to flow rates and other rate dependent param- 
eters entering the mass balances. Likewise, in diffusional systems the local 
state is determined by mass balance relations associated with mass transfer 
processes, in the following subsections the mass balance relations appro- 
priate to various systems are summarized. 

4. 2. 3.1 Expansion of isolated Systems 

in the expansion of a fixed mass, closed, adiabatic system it is usually 
appropriate to trace its state history as a function of static pressure. If 
the process is reversible, the entropy is constant and the local state is not 
a function of the time history of the expansion. Such systems satisfy the 
simple mass balance constraint 

a . = constant (65) 

i 

This equation implies either total equilibrium or a mixed equilibrium- frozen 
chemical process. If, however, finite reaction rates are important, the path 
ceases to be reversible, entropy rises, and the time history of the expansion 
must be considered. If the pressure is a known function of time, the expan- 
sion can be treated as 


state = f(a^,s,P) (66) 

where the state includes such terms as da^/dO and ds/d0. The rate of change 
of the "elemental" composition is obtained from Eqs . (57) and (58), whereas 

the rate of increase in entropy is given by 


ds 

dO 



In K 


E 



p j 


m P 771 


(67) 


This derivative is well behaved, even as equilibrium is approached and may be 
evaluated explicitly if desired. The derivatives, however, must be 
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treated implicitly if any hope for near equilibrium solutions is to be main- 
tained. Once a particular formulation is adopted the techniques suggested in 
the nonequilibrium presentation (Section 4.2.2) can be introduced in order to 
assure consistent solution validity. Because of the simplicity of the mass 
balance relation for this process, it is practical to include the kinetic 
mass balances directly with the iterative solution of the chemical state. 

This state calculation includes the relations previously presented together 
with the entropy constraint, namely 



3 


= P^s 


( 68 ) 


4. 2.3.2 State Calculations for Open Systems 

The evaluation of the state in a general open system involving diffusive 
and convective mass and energy fluxes is most generally performed as a sub- 
ordinate solution. For example, in the boundary- layer solutions of current 
interest, state solutions are required at several interacting locations, in 
this application, state solutions are required based on assigned "elemental" 
mass fractions (or cu) , enthalpy and pressure, i.e. 


state = f(Gu,h,P) 


(69) 


This solution provides to the boundary-layer solution the detailed state in- 
cluding thermodynamic, transport and radiative properties as well as the pro- 
duction rates da^/dt. The last term must be included with the general mass 
balance relations of the boundary -layer program. 

The specific relations used to achieve the state solution are the equi- 
librium equations (Eq. (48)), the mass balance relations (Eq. (54)), the pres- 
sure constraint (Eq. (51)) and an enthalpy constraint 



(70) 


which involve no greater complexity than the conventional isolated system 
equilibrium solution. The coupling between this solution and the boundary- 
layer solution requires not only the evaluation of the production rates, 
da^/dt but also the rate of change of these rates with respect to the inde- 
pendent parameters on the right-hand- side of Eq. (69) . The rates are deter- 
mined with Eqs . (58) and the derivatives by use of relations developed in 
Part V of this series of reports. 


Again, the problem of stability threatens when equilibrium approaches. 
However, by following the approach previously presented, the kinetic terms 
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can be treated separately While all other terms of the mass-balance equations 
are being collected. These equations can then be rearranged and combined with 
the kinetic relations in such a way that controlling reactions again affect 
only one equation at each location. Because of the overall implicit character 
of the boundary- layer solution, this procedure will, on convergence, yield 
valid stable solutions. it has been found, however, that the introduction of 
equilibrium type relations into the set of boundary-layer equations can de- 
stroy the linearity of the system. Therefore the approach of equilibrium by 
the kinetic equations included in the boundary- 1 layer mass balances would 
probably delay convergence and necessitate the inclusion of certain itera- 
tion constraints. 


4. 2. 3. 3 Surface State Solutions 

A more complex set of mass balance relations are introduced when surface 
state .solutions are sought. Coupling between boundary layer, internal conduc- 
tion, and surface mechanical removal solutions may be involved in these rela- 
tions. in effect all the other mass balance solutions become subordinate to 
this solution. Two types of boundary-layer representations have been devel- 
oped, a transfer coefficient correlation of mass transfer using the Z-potential 
and the expression of the fluxes at the wall in a form compatible with the 
boundary layer nodal solution procedure. The governing equations, including 
heterogeneous reaction kinetics, are presented in Part V of this series of 
reports . 

One of the most elusive aspects of surface-state solutions is the adequate 
specification of the mechanical-chemical surface constraint. The present for- 
mulation is based on the following set of constraints for condensed phase 
species . 

p = 0 if T < T (71) 

1 i 

and 



Xn K 

P X 


(72) 


with the equality applying to one species with Tp^ s T. The first equation 
implies that a particular condensed species cannot leave the surface until 
the surface temperature is at or above that species fail or flow temperature, 
Tr. . The second relation states that all present condensed species are in 

* j i 

equilibrium with the base species. The inequality applies to non-present 
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condensed species and prohibits a super-saturated vapor state. This is equiva- 
lent to saying that at 100° C 


P H 2 0 — 1 atm 


(73) 


The requirement that one species be at or below its fail temperature estab- 
lishes the structural limitation of the surface. A typical result might 
show a surface at 2500°K with the equilibrium holding for SiC* and SiO*, but 

if SiO? has been assigned a fail temperature of, say, 2300°K, p_ . will be 
z Si0 o 


positive indicating liquid removal of SiO*. The SiC*, with a fail temperature 
greater than 2500° K, represents the surface constraining species. 


4.2.4 Oblique Shock Relations 

In the case of an oblique shock wave, constraints of Eqs. (51), (68), 

and (70) do not apply and are replaced by the conventional equations for 
conservation of energy, mass and momentum across the oblique shock. Presum- 
ing knowledge of the upstream conditions, these relations yield as unknowns 
only h, p, and P downstream of the shock. These relations can be further 
reduced to 


C— ■ , H m (p-u.coseo 2 

+ <»!»!«*•<>!>* = P 1 + ^ 

j 

Y p i H j + i <pi u i cose i > 3 = h i + uic r 9i) ] 


(74) 


(75) 


where the non-subscripted variables are downstream of the shock. The first 
of this pair of equations replaces the more conventional pressure constraint 
and the latter the enthalpy constraint. 


4.2.5 Summary 

In this subsection an attempt has been made to formalize the basic rela- 
tions so as to simplify the generation of an orderly solution. Unfortunately 
dealing with nonlinear equations such as these is never straightforward and 
is subject to many pitfalls, in the next subsection, the procedures as adopted 
in the current solution technique are described. 


4.3 NUMERICAL SOLUTION PROCEDURE 

The solution to a set of simultaneous nonlinear algebraic equations can 
be either trivially simple or agonizingly difficult, depending on the linearity 
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of the system and the depth of coupling existing between the equations. None 
of the problems formulated in Section 4.2 fall into the first class and some 
fall into the latter. The basic formulation adopted is relatively conventional 
and will be described first, followed by some discussion of the pitfalls that 
can be encountered and devices adopted to circumvent them. 

4.3.1 Basic Formulation 

The most direct method of solving a set of nonlinear algebraic equations 
is the Newton-Raphson procedure. Its application is straightforward in con- 
cept but in reality many choices occur during the formulation of a specific 
problem, choices which can affect the success or failure of' a specific solu- 
tion. The method itself is the extension of Newton’s iterative method to 
multi-dimensional problems. Errors are evaluated for each of the equations 
based on a set of trial values for the unknown independent variables. The 
rates of change of these errors with respect to these independent variables 
are analytically determined, also based on the trial values. In the formula- 
tion which has been followed, in p., p. , in T, and in (Pftft are taken as 

J x 

the set of independent variables, but corrections are often in terms of p j , 

1/T and P 771 , In some systems this choice yields linear mass balance equa- 
tions which if once satisfied will never deviate. 

4.3.2 Solution Convergence 

In general, the convergence of the set of equations appropriate to a 
particular problem depends on a number of factors in addition to the formula- 
tion of the derivatives. in addition to the selection of correction coordi- 
nates, initial estimates and correction restraints are major factors. 

4.3. 2.1 Correction Restraints 

In highly nonlinear application of the Newton-Raphson technique, a variety 
of constraints with regard to independent variable corrections are necessary. 
These constraints all manifest themselves in a damping factor which limits 
the extent which the solution is advanced down the correction vector. When 
corrections exceed the constraint limits, a damping factor is introduced which 
is applied uniformly to all variables. 

4.3.2. 2 Initial Guesses 

It is obvious that a good first guess can save time in any iterative solu- 
tion. In the present formulation these guesses are generally based on previous 
solutions and only the initial stagnation or shock solution does not have the 
benefit of prior solution. This solution is readily obtained from practically 
any first guess, since the stagnation state is usually at relatively elevated 
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temperatures and has a fixed "elemental" composition, in the subroutine ver- 
sion of the chemical state program used in conjunction with the boundary-layer 
procedure, first guesses are generally based on solutions at the same boundary 
layer transverse location stored during prior iterations in the boundary-layer 
program or from solutions at the preceding axial station. 

Because of the introduction of new species by the wall material it is 
necessary to initialize their compositions when the corresponding elements 
appear in the state solutions. Likewise if a species disappears, e.g., as 
the edge of the boundary layer is approached or because the sequence of bound- 
ary-layer iterations results in the termination of surface mass addition, it 
is necessary to zero the species in a fashion that will not result in a singu- 
lar solution for the rest of the equations. 

Bookkeeping becomes a major factor in the state programs if efficient 
and stable repetitive utilization is to be made of the routines. This book- 
keeping establishes optimum first guesses, determines which atomic elements 
are present and zeros or initializes the appropriate molecular species. 

4.4 CHEMICAL STATE PROGRAMS 

To treat the solution of chemical state problems, two equilibrium pro- 
grams are currently in operation, namely, the equilibrium surface thermo- 
chemistry (EST) program and the Aerotherm chemical equilibrium (ACE) program. 
For nonequilibrium systems a special version of the ACE program is utilized 
with the KINET subroutine. Ultimately, this latter combination will be gen- 
eralized into the general nonequilibrium ablation thermochemistry (GNAT) pro- 
gram. In this section, the current capabilities of these routines is sum- 
marized, together with brief descriptions of the input required and output 
obtained for operation under various program options. 

4.4.1 The Aerotherm Chemical Equilibrium (ACE) Program 

The ACE program is the more recent and more general of the two equilib- 
rium programs . The program operates either as a separate routine or as a 
subroutine to the boundary layer integral matrix procedure (BLIMP) . Some 
options apply in both cases, but certain bookkeeping aspects are modified 
in order to streamline the subroutine version. The principal options are as 
follows : 

1. Oblique shock relations 

2. Assigned enthalpy, pressure and elemental composition 

3. Assigned entropy, pressure and elemental composition 

4. Assigned temperature, pressure and elemental composition 
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5. Surface balances with assigned temperature or assigned 
surface equilibrium 


The equations actually involved in these calculations have been summarized 
in the preceding subsections and developed in Part V of the current report 
series. All of the options are formulated for a general chemical system. 
Consideration of any molecular, ionic or atomic species requires only the 
inclusion of the basic thermochemical data appropriate for that species. 

These basic data are obtained, for example, frqm the JANAF Thermochemical 
Data Tables and include entropy, heat of formation and specific heat curve- 
fit data. 

Of the options listed, only the first and last require further definition 
of the requisite input. The oblique shock option accepts the upstream veloc- 
ity, density and static enthalpy and the shock angle as basic input along with 
the elemental composition. The output includes the state of the gases down- 
stream of the shock and the isentropic stagnation state. 

The surface mass balance options require as a minimum input the normalized 
pyrolysis gas and char recession rates as well as the elemental composition of 
these components and the pressure. If surface equilibrium is to be suppressed, 
temperature must also be assigned. Two forms of surface mass balances are 
included in the ACE program. For coupling with the BLIMP program, special 
linearized flux relations are developed by a truncated Taylor series expan- 
sion about the current trial wall fluxes and wall state (see Section 2.3.1) . 

For use with transfer coefficients, the program requires the specification 
of edge composition and, further, if the unequal diffusion model is to be 
used, the diffusion factors, F i# must be specified unless the logarithmic 
proportionality of these factors to molecular weight is utilized (see Section 
2.1.2). When liquid-layer removal is contemplated, it is necessary to specify 
the maximum temperature at which each condensed species can structurally sup- 
port the surface. The output from the surface mass balance options include 
the total definition of the surface state including temperature, condensed— 
material removal rate, and the condensed species which structurally maintains 
the surface. Output from this option is also obtained on cards suitable for 
direct use with the CMA program if transfer coefficient mass balances are 
performed. 

For most options a rather complete state of the system is generated which 
includes compositions, thermodynamic and transport properties, and major prop- 
erty and composition derivatives. It is these derivatives which permit the 
analytic treatment of the complex boundary- layer equations in the BLIMP program. 
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4.4.2 The Equilibrium Surface Thermochemistry (EST) Program 

The EST program has the advantage of seniority and experience but in- 
cludes no options excluded from the ACE program. Thus the advantages of EST 
involve better diagnostic output, more formalized input and output, and, on 
some occasions, greater convergence stability. EST is designed primarily for 
the surface mass balance options but also can be used for assigned tempera- 
ture, elemental composition and pressure solutions. It does not include the 
condensed phase removal model, and is limited to the transfer coefficient 
mass balance relations. A typical output (from Ref. 20) is shown in Fig. 8. 

4.4.3 The ACE-KINET Program 

For nonequilibrium solutions, the KINET subroutine adds to the surface 
mass balance option of ACE the ability to treat specific heterogeneous or 
homogeneous surface-catalyzed reactions. Special KINET routines are prepared 
for specific systems and include a predetermined set of kinetically controlled 
reactions. The routine currently in use treats the heterogeneous oxidation 
of graphite by CC^# 0 ^, ^O, and the surface catalyzed water gas shift reac- 
tion in the H, C, 0, N system. A typical result for graphite phenolic abla- 
tion (Fig. 9) shows the significant low-temperature effect of the kinetic 
relations. The basic data required for specifying the kinetic rates of each 
reaction are activation energy, pre-exponential factor, and reaction order. 

The output is identical to the ACE program but does not include nonequilibrium 
state derivatives. 

4.5 SUMMARY AND CURRENT STATUS 

In the preceding subsections, an approach for determining the equilibrium 
or nonequilibrium chemical state has been summarized for a number of open and 
closed thermodynamic systems. An effort has been made to provide a relatively 
general approach to the problems associated with such solutions and to indi- 
cate means of circumventing them. A brief discussion of the mechanics of 
the solution served to introduce the program and subprograms involved in the 
computer analysis. Some of these routines are quite general in their present 
formulation, others are directed toward specific systems. 

Currently all equilibrium aspects of the program are fully operational 
for general chemical systems. This includes the various closed and open sys- 
tem options, the shock wave relations, the surface coupled boundary layer mass 
balances, bookkeeping involved with treating appearing and disappearing atomic 
elements, and the property and property derivative calculations. The KINET 
routine currently treats only the heterogeneous reactions associated with 
graphite oxidation. The generalization of this routine following the approach 
presented in this section is a major recommendation of this report. 
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ABLATING SURFACE THERMO-CHEMICAL EQUILIBRIA 


PAGE l 


PREPARATION OF 
GENERATOR 


TABLE FOR GRAPHITE PHENOLIC ABLATION IN ARC PLASMA 

7/6/64 


RELATIVF ELEMENTAL COMPO S I T I ONS » ATOMIC WTS/UNIT MASS 

AT. NO* ELEMENT ATOMIC WT PYRO.GAS CHAR LAYER 

1 HYDROGEN I.OOBOC 0.1070182 *0. 

2 HELIUM 4.00300 -0. -0. 

6 C ARBTIN 12.01100 0.0505157 0.0832570 

7 NITROGEN L4. 00800 -0. -0. 

8 OXYGEN 16.00000 0.0178364 -0. 


B/L EDGE GAS 
- 0 . 

0.0570572 

- 0 . 

0.0441962 

0.0095312 


PRESSURE,ATM 9.00000E 00 SURF • TEMP f K 3559.68 ENTH,CAL/GM 1.78067E 03 
MOOT P.G./CM 5 • OOOOOE-02 MOL. WEIGHT 13.2584 HWU+M/CM) 2.22583E 03 
MOOT CHAR/CM 2.00000E-01 SPEC. HEAT 0.51900 SURFACE C* 22 

- - - - MOLE FRACTION - SPECIE - - - - 

0.0005308 C 0.0000178 CH 0.0161998 CHN 0.0000000 CHNO 

0.0000035 CHO 0.0000004 CH2 0.0000000 CH20 0.0000004 CH3 

0.0000000 CH4 0.0045728 CN 0.0005663 C2 0.0064974 C2H 

0.0009285 C2H2 0.0012684 C2N2 0.0000001 C2H3 0.0000000 C2H4 

0.0000000 C2H40 0.0000000 C2H6 0.0019618 C3 0.0060477 C 3H 

0.0000228 C3H2 0.0000002 C3H3 0.0000000 C3H4A 0.0000000 C3H5 

0.0000000 C 302 0.0000331 C4 0.0063986 C4H 0.0000481 C4H2 

0.0000000 C4H3 0.0000000 C4H4A 0.0001094 C4N2 0.0000493 C5 

0.0000174 C5H 0.0000003 C5H2 0.0000000 C5H3 0.0000010 C6 

0.0003574 C6H 0.0000016 C6H2 0. 0000000 C6H3 0.0000000 C6H6 

0.0000007 C 7 0.0000077 C7H 0.0000000 C7H2 0.0000000 C8 

0.0000000 C 8H 0.0000000 C8H2 0.0000000 C9 0.0000001 C9H 

0.0000000 C9H2 0.0000000 CIO 0.0000000 ClOH 0.0000000 C10H2 

0.0127292 H 0.0000041 HN 0.0000000 HNO 0.0000000 HN02 

0.0000000 HNO 3 0.0000000 HO 0.0000000 H02 l 0.0000457 N 

0.0000000 NO 0.0000000 N02 0.0000000 N20 0. 0000001 0 

0.0000000 02 0.0000001 C02 0.1105511 CO 0.6051916 HE 

0.0032347 H2 0.0000000 H20 0.2225999 N2 


1 

| 


i 

i 

l 


Figure 8. Typical Output From Equilibrium Surface Thermochemistry 
(EST) Program 
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Thermo chemical Response 



The section has discussed in rather general fashion the treatment of 
general chemical systems. The ultimate program which should evolve from this 
analysis will be a General Nonequilibrium Ablation Thermochemistry (GNAT) pro- 
gram designed for treating the problems associated with equilibrium and non- 
equilibrium at and above ablating surfaces. 


SECTION 5 

COUPLED BOUNDARY LAYER AND CHARRING MATERIAL 
COMPUTATIONAL PROCEDURES 


5.1 INTRODUCTION 

The charring material response computation described in Section 3.3 above 
can, of course, be operated independently of any boundary— layer calculations 
if the temperature and recession rate are specified as the known boundary con- 
dition. In fact, the CMA program has many uses when operated in this way, 
particularly as a tool for extracting material properties (thermal conductivity, 
specific heat, and decomposition kinetic parameters) from measured test data. 

To make a general predictive tool, however, the CMA program must be cou- 
pled to some boundary- layer calculation. For best accuracy, this calculation 
should be complete in all relevant details. The BLIMP program described in 
Section 2.3 constitutes such a boundary- layer procedure. The coupled version 
of BLIMP plus CMA, denoted the CABLE program, thus provides a complete charring 
ablator analysis procedure. This coupled procedure is discussed in Section 
5.3 below. 

For greater speed, the CMA program can be coupled, instead, to the ACE or 
EST programs described in Section 4.4 which approximate the boundary layer by 
convective transfer coefficients but still retain the essential chemical fea- 
tures of the ablation events including the effects of unequal diffusion coef- 
ficients for all species. Section 5.2 below describes the approach utilized 
for obtaining solutions with these programs. 

5.2 CHARRING MATERIAL ABLATION PROGRAM COUPLED TO FILM-COEFFICIENT 

BOUNDARY- LAYER CHEMISTRY PROGRAM: EST OR ACE 

5.2.1 General Problem Description 

Section 3.2.3 above describes how the CMA in-depth solution routine can 
be coupled to a surface energy balance procedure to provide the heated surface 
boundary condition. The surface energy balance was given in normalized form 
in Eq. (30) . It can also be expressed as 
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q rad 
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( 76 ) 


where 


q rad ” q rad 
in out 


(pv) w = ra g + * 


* C - Zv 


(77) 


(78) 


The m g and q cond - f ( T vP are delivered by the in-depth solution. Other 
dependencies of interest are 


h g - V T vP 


h c = W 


^rad 

out 


q rad^ T w^ 

out 


For the other terms 


-q a 


q rad' 

in 


2 X 


h, = 


functions of boundary- layer- edge enthalpy, 
pressure, local boundary- layer solution, 
laws for conservation of chemical elements, 
chemical equilibria and/or kinetic rela- 
tions, upstream events, m and m 

c g 

From the standpoint of the surface energy balance solution the desired 
relationship may be summarized as 


-q a 


q rad' h w 
in 


- E 


m h. 


_ function? of m c 

for given P 


and rft 


(79) 


Equation (76) requires an iterative solution in which T and m are the 
, t t w c 

primary variables of interest, one of them regarded as dependent and the 

other as independent. It is most convenient to obtain the relations of Eq. 
(79) outside of the CMA in-depth solution and to provide the resulting in- 
formation as tables which may be stored and referred to as needed. These 
tables give -q-^, q ra(J in , h w , ^ iti^h^, and T w as functions of m,, 

Ag, and another variable (essentially time but including all time-dependent 
aspects such as pressure) . 
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The energy balance solution procedure proceeds as follows. An initial 
i of the char consumption rate, m,, is obtained in some manner. With 
m and the m supplied by the in-depth solution, the quantities 


~ q a w # q rad in' 


h w - Z 


and T. 


are obtained by table look up in the 


tables provided by the outside surface solution routine. The quantities 
b , h and q^^ Qut are then formulated using the T w so obtained. The 
surface energy balance (Eq. (76)) is then computed, the q cond as a function 
of T having beeh provided by the in-depth solution. In general, however, 
the sum of the terms will not equal zero but some error. An iteration pro- 
cedure is then used to select successively better estimates of m c which drive 
the error to zero. Experience shows that Newton's procedure, in which the 
derivative of the error with respect to m c is used to compute the next guess 
for m c , gives good results. 

The following sections describe how the required tables are generated 
using the film-coefficient model of the boundary layer. First, the film 
coefficient expressions for -q a are presented, followed by a discussion of 
the chemistry-and-mass-balance solution required for constructing the surface 
tables . 


5.2.2 General Requirements for Energy Flux 

Film coefficient correlations of boundary-layer heat-transfer analysis 
and data have been used for years to predict mass and energy transfer rates 
at a surface. If no chemical reactions are involved, the adaptation of the 
film-coefficient model to the mass -transfer problem is rather simple and 
straightforward. If equilibrium chemical reactions are involved, either in 
the boundary layer or at the surface, the ultimate form of the film-coeffi- 
cient adaptation remains relatively simple for the special case of equal mass 
diffusion coefficients for all species and unity Lewis number. In this case, 
the diffusive energy flux -q is given by the familiar expression 


-*a = Pe u e C H (H r " V 


(80) 


where C H is the Stanton number for heat transfer. This expression has been 
generalized to include the effects of and unequal mass diffusion 

coefficients as described in Part II of the present series of reports. in 
that case the expression for -q^ becomes more complex and involves more in- 
formation about the wall state than simply h w - The extra information auto- 
matically is computed in the process of computing the wall state in any case, 
as noted in the next section, so no extra computation is involved for the 
more complex case. 
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5.2.3 Solution Procedure as Related to Tabular Formulation 


Construction of the tabular boundary condition to the CMA program basi- 
cally requires solution for the state of the gas adjacent to the heated wall, 
given the normalized char rate s m c /p e u' e C^ and the char chemical compo- 
sition, the normalized gas injection rate s ^ l g y ^^e U eSl an< ^ 9 as c ^ ein ^“ 
cal composition, and the boundary- layer edge state.. Such a solution provides 

h and -q a needed for the surface energy balance. The required solution 
w w 

satisfies the necessary mass balances (with mass diffusive fluxes given by 
film-coefficients analogous to those for the diffusive energy fluxes) and 
chemical relations, but of course does not attempt to satisfy an energy bal- 
ance, since this last step is to be performed by the CMA program using the 
tabulated results of the chemistry-plus-mass-balance solution. 


Several computer programs are available for the determination of the 
wall state tables. The Equilibrium Surface Thermochemistry program (EST) 
is one such program for the special case of chemical equilibrium at the wall. 
The Aerotherm Chemical Equilibrium program (ACE) is a later version of the 
EST program which allows the surface run-off of materials formed at the 
surface above their fail temperatures. Hence ACE computes ]Tm r ^hg and lumps 
it with -q_^ for tlie later solution of the surface energy balance. Utiliza- 
tion of the KINET subroutine together with the ACE program allows heterogeneous 
kinetic control of several simultaneous reactions with the surface material 
as well as surface-catalyzed homogeneous reactions. These routines have been 
discussed in Section 4.4. 


The chemistry-plus-mass-balance solutions provided by the EST or ace pro- 
gram supply as output tables of T^, h^ and other quantities needed to com- 
pute -q , all as functions of B', B' # and pressure. During each time step 

, c *w c y 

m the course of the in-depth response solution, the CMA program develops an 

expression q CQnd = ( T w ) » substitutes this into the surface energy balance 

(Eq. (76)), and then searches among the surface tables for a B 1 which yields 

an energy balance, thus defining a new value for T . Then the solution pro- 

w 

cedure is ready for the next time step. 


It may be noted that the tabular approach to the surface chemistry solu- 
tion is suggested by economy. Without such tables each iteration in the search 
for a surface energy balance would require a new surface chemistry solution, 
generally in the near neighborhood of many such previous solutions, in almost 
all problems, the tabular approach involves fewer surface computations. Fur- 
thermore, tables once generated are often useable for many different problems, 
yielding even greater economy. Finally, without the tabular approach, the 
occasional nonconvergent surface chemistry solution would stop the entire in- 
depth solution process, with the tabular approach, such solutions are auto- 
matically weeded out of the tables without damage to the subsequent in-depth 
solution . 
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5.2.4 Experience with Film-Coefficient Coupling 

The in-depth CMA solution routine coupled to the film-coefficient models 
provided by the EST, ACE, and ACE/KINET programs has been tested on a fairly 
wide variety of materials and a brief account of these tests may provide some 
useful orientation for the reader. (It may be noted that it is sometimes 
difficult to establish whether any discrepancies between predictions and data 
are due to errors in in-depth properties or models or errors in the surface 
treatment. Usually a careful testing plan first establishes whether or not 
the in-depth aspects are being treated correctly by comparing in-depth pre- 
dictions to in-depth data, such as thermocouple response, where the predic- 
tions are made with assigned surface temperature and recession rate to match 
observed surface data. If the in-depth model can be verified in this manner, 
calculations may then be done with the general thermochemical boundary con- 
ditions. Success in predicting surface temperature and recession under these 
circumstances constitutes the "good results" referenced to below.) 

By suppressing pyrolysis effects, the program has been used for numerous 
transient ablation problems featuring noncharring refractories. Examples 
here have included alumina, boron nitride, tungsten, and graphite and have 
covered problems with liquid layer runoff and kinetic control with generally 
excellent results . 

With regard to charring materials, the program has been run very exten- 
sively for graphite-phenolic and carbon-phenolic with good results. Success 
with nylon-phenolic has been mixed since this material often suffers from 
mechanical ablation effects not included in the program model; the same re- 
marks apply to asbestos-phenolic. Materials with substantial silica content 
have been frequently predicted, sometimes with good success, but other times 
with poor results since materials with large amounts of silica occasionally 
display such physical events as thick liquid layer runoff and subsurface char 
reactions (for example, silica-carbon reactions) not accounted for by the 
program . 

In conclusion, the coupled computation procedure constituted by the CMA 
program plus some film-coefficient based chemistry solution (EST, ACE, or 
ACE/KINET) has been applied to a wide range of materials of technical interest 
with excellent to poor correlation depending on the particular material and 
boundary conditions. Any discrepancies between predictions and data have been 
clearly attributable to effects not considered by the program or occasionally 
to ill-judged boundary conditions or material properties. The program appears 
to be fully checked out and operational for the physical and chemical models 
currently employed. 
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5.3 COUPLED ABLATOR/BOUNDARY LAYER/ENVIRONMENT (CABLE) PROGRAM 


5.3.1 Introduction 

The coupled ablator/boundary layer/environment (CABLE) program is a com- 
putational procedure which couples the transient response of a charring heat 
shield material to a chemically reacting laminar boundary layer appropriate 
to superorbital reentry. A coupled approach is necessary since the material 
response affects the structure of the boundary layer, and the boundary layer 
determines the energy and mass fluxes at the surface which in turn control 
the heat shield response. The CABLE program incorporates subroutine versions 
of the BLIMP program (described in Section 2.3) and the CMA program (described 
in Section 3.3) . The features of the CABLE program are summarized in Section 
5.3.2, the mechanics of coupling are discussed in Section 5.3.3 and the cou- 
pling procedure is demonstrated further by a sample problem presented in Sec- 
tion 5.3.4. 

5.3.2 Characteristics of the CABLE Program 

All of the features of the BLIMP and CMA programs pertinent to the cou- 
pled problem are retained in the CABLE program. In that the characteristics 
of these subprograms have been described in some detail in preceding sections, 
the models employed in the CABLE program are presented summarily in Table I. 
The operational status of the various aspects of the computational procedure 
are also summarized therein.. It can be seen that many considerations are 
fully operational, including all aspects of the in-depth response of the ablat- 
ing surface material and the nonablating backup material. Certain aspects of 
the boundary-layer solution cannot be considered fully operational until such 
time that the procedure is checked out for the wide variety of materials, 
environments and flight conditions for which it is presumably applicable. 

Some aspects of the ultimate boundary- layer program have not been fully imple- 
mented at this time. Those areas where additional effort is recommended are 
discussed in Section 6. 

5.3.3 Coupling Procedure 

In the present approach, a series of one-dimensional transient charring 
ablation solutions are directly coupled to time-varying but quasi-steady two- 
dimensional boundary layers as shown in the following sketch. From a study 
of the numerical equations associated with the boundary- layer and charring- 
ablation solution procedures, it is seen that the charring-ablation solution 
at station l and time 0 is dependent upon the boundary- layer solutions at 
previous stations l-l and i-2 through the three-point finite difference 
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IV IN-DEPTH RESPONSE OF EXPOSED 

(ABLATING) MATERIAL (Continued) 
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relations for the streamwise derivatives, is dependent upon the 0 - A0 solu- 
tion of the internal response of the ablation material at & through the 
charring ablator finite-difference relations, and is implicitly dependent 
upon itself and the current boundary-layer solution. 

Several approaches for coupling the boundary layer and charring ablation 
solutions were considered in the present study. The method which was finally 
adopted was selected on several bases: it makes use of options of existing 

programs which are well exercised and known to perform well, it avoids extra- 
polation of surface boundary conditions, and it avoids repeated (iterative) 
solution of the boundary layer and transient charring ablation response. 

The various other methods which were considered were inferior in one or more 
of these considerations. Furthermore, storage requirements and computational 
time are improved relative to most if not all of the other methods considered. 

In the procedure which has been adopted, the transient charring ablation 
solution is effectively the controlling program. The charring ablation solu- 
tion at a given station proceeds noniteratively , calling the boundary-layer 
procedure as needed to supply the surface boundary condition. The complete 
transient history at each axial station is performed prior to advancing to the 
next axial station. This is accomplished by performing sets of nonsimilar 
boundary-layer solutions at the current station for a discrete array of times 
(0) , normalized pyrolysis-gas mass flow rates (m*) , and normalized char mass 

y 
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flow rates (m*) (or surface temperatures, T^, when m* = 0) which bracket the 
current values for thes'e parameters. Calculations for intermediate times and 
intermediate values of m* and m* (or T^) are then performed by interpola- 
tion as they are needed for the charring ablation solution. It is significant 
that only those members of the m*, m* array needed to contain current 
values are considered, and that at any instant, these are required for only a 
pair of times. 


The procedure is demonstrated by the example illustrated in the follow- 
ing sketch which is a planar representation of three-dimensional (m*, m*, 0 ) 
space. 


m* 

g 



m* 

c 


m* 


c 


The lines of arrows shown in the sketches are the projections in the time 
planes 0=0^ and Q = 0^ of a hypothetical history of m* and m* in a 
coupled boundary layer and charring ablation solution between times 0 ^ and 
where 6.^ > 0^. The solution at time 0^ is indicated by asterisks, 
whereas the solution at time $ 2 i s indicated by circles. The times 0^ and 
and the grid values for m* = 0,1,2,... and m* = 0,1,2,... are preselected 
values for these parameters at which parametric boundary- layer solutions are 
conducted if and when needed. Based on the point (*) at time 0^, boundary- 
layer solutions are generated for the m*, m* points 1,1; 1,2? 2,1; and 2,2 
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at times 0 ^ and 0 2 .* Charring ablation solutions can be obtained for times 

0 1 <;0<:0 2 by linear interpolation as long as m* and m* stay within 

these values . Suppose that the course of the calculation between times 0^ 

and 0 2 are as indicated in the sketch. Then, additional solutions at m*, 

m* of 1,3 and 2,3; then 3,2 and 3,3; and finally 2,4 and 3,4 would be re- 
c 

quired, each at both times. When time 0 ^ (point @) is approached, the BLIMP 

program is called upon for solution at time 0 ^ for the exact values of m* 

and m* . This boundary- layer solution is printed out and that information 
c 

needed for future reference (at downstream stations) is saved on tape. Solu- 
tions are then performed for time 0 for the current bracketing values of 

m*, m* (in the present example, values of m*, m* of 2,3; 3,3; 2,4; and 
g c 9 c . _ 

3,4). These boundary-layer solutions at time 0^ are placed over those for 

0i a tape flip-flop since the latter are no longer needed. The charring 

ablation solution next proceeds from time 0 2 to time 0 3# calling the BLIMP 

program only in the event that this range of m*, m* is exceeded. 

By the use of this procedure, storage requirements are nominal. In the 
first’ place, the charring ablation solutions are noniterative and the complete 
transient solution at a station is accomplished and the results printed out 
before advancing to the next station. Thus no historic information relative 
to the charring ablation solution has to be stored. With regard to the bound- 
ary layer, only two times with four m*, m* combinations at each of these 
times need be considered at the same time. The only quantities in the boundary 
layer which need to be dimensioned for the full time array are three input 
quantities of time, total pressure, and total enthalpy. Edge conditions are 
computed around the body at the time of the stagnation-point calculation since 
the necessary integrations are performed by curve fitting. This necessitates 
that streamwise dimension, static pressure, edge velocity, edge density, edge 
viscosity, edge temperature, body curvature parameter ( r Q ) * transformed stream- 
wise dimension (T) * pressure gradient parameter (g) , and the flux normalizing 
parameter (a*) be dimensioned for the number of streamwise positions (but not 
for time) . About 300 numbers must be stored during the flip-flop operation 
associated with the two times which are being considered simultaneously, whereas 
about 500 numbers must be stored on tape to reenter the boundary layer at the 


*This simple example applies only after char recession has commenced, that is, 
m* > 0. In the early portion of a trajectory when m* = 0, m c is replaced 
by T as an independent parameter. Furthermore, it is then necessary for 
each W 0 and mi of current interest to compute the T w above which char 
recession would occur. This is illustrated in the sample problem presented 
in Section 5.3.4. 


81 



same time but at the neSt downstream station (used for first guesses and for 
calculation of nonsimilar terms) . Thus, both permanent machine storage re- 
quirements and tape storage requirements are not excessive as a consequence 
of coupling. 

This coupling approach has the important feature that the CMA program 
operates very nearly as it does when used in conjunction with the ACE program 
(see Section 5.2). In the CMA/ACE approach, complete surface tables are com- 
puted a priori and independently with the ACE or EST program and these are 
available to the CMA solution. In the coupled approach, these surface tables 
are initialized with the word VOID. When the CMA program encounters this 
word, the BLIMP program is called to supply the requisite information for 
that 9, m* and m* (or T w ) . It is thus significant that the CMA/EST and 
CMA/ACE approaches have been used extensively and very successfully for a 
wide variety of materials and environments. Likewise, the boundary-layer 
calculations are performed with assigned m* and m* or assigned m* and 
V to ^ ether with the requirement of surface equilibrium (with possible spe- 
cified rate-controlled surface reactions) , options of the BLIMP program which 
also have been exercised extensively with success. Furthermore, this replace- 
ment of the wall mass and energy balances by these simple assignment state- 
ments adds stability to the boundary- layer solution. 

5.3.4 Coupled Solution for Apollo SA 202 Trajectory 

As a demonstration of the coupling procedure, a trace of m* versus 
T w at the stagnation point of the Apollo heat shield during the^irst 76 
seconds of the Apollo SA 202 reentry trajectory is presented in Figure 10. 

In this problem the time-table entries were selected to be 4310, 434Q, 4375 
and 4400 seconds; m* entries as 0, 0.1, 0.2, 0.3; T w entries as 500, 1000, 
1500, 2000, and 2500 R; and m* as 10" 5 , 10 -3 and 10~ 3 , Boundary- layer solu- 
tions were performed at combinations of these independent parameters as they 
were required and are numbered in the sequence in which they were performed. 

The first step in the coupled solution was to initialize the charring 
ablation solution at the assumed initial temperature of 530°R. This was 
followed by a boundary-layer solution at this initial nonablating condition 
(0^ = 4310 sec, T w = 530°R, m* = 0) . This is identified as Solution 1 in 
Figure 10. The next step was to find the wall temperatures at which ablation 
would start for the initial time of 4310 seconds and the second entry in the 
time table, 4348 seconds, each for the first two entries in the m* table, 
namely 0 and 0.1 (Solutions 2 through 5). This was accomplished by computing 
the surface temperatures required to maintain surface equilibrium for these 
boundary-layer edge conditions (i.e., times), these normalized pyrolysis gas 
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U3* 


Nonablating 


Ablating 






Wall Temperature °R 

c 

Figure 10. Demonstration of Coupling Procedure for Apollo 
Stagnation Point During SA 202 Trajectory 
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flow rates, and a very sjjiall normalized surface recession rate (m* of 10 -5 
was used in these boundary-layer calculations) . it can be seen that these 
n ablation temperatures" are all higher than the first two entries in the T 

w 

table of 500 and 1000 R. Therefore, T w is the appropriate independent param- 
eter (rather than m*) during this portion of the trajectory. Boundary-layer 

solutions were then obtained at the eight corners of the 0, m*, T cube 

g w 

(Solutions 6 through 13) . At this point the transient charring-ablation solu- 
tion was able to commence, interpolating between the bracketing values of 0, 

m* and T . 
g w 

The transient charring-ablation solution then proceeded to 0^ = 4348 
seconds, the time steps being determined by various controls' built into the 
implicit finite-difference procedure, A charring-ablation solution was per- 
formed at precisely 4348 seconds, and this was followed by a boundary-layer 

solution at that time, wall temperature, m* and m* (Solution 14) . The next 

5 c 

order of business was to obtain the ablation temperatures (Solutions 15 and 
16) at the next entry in the time table, 0^ = 4375 seconds, for the two cur- 
rent m* of 0 and 0.1, and to obtain boundary- layer solutions at this 
new time for the two current T^, and m* (Solutions 17 through 20) , 

The transient charring-ablation solution was then able to recommence and 
continue until a surface temperature of 1000°R was attained. It was then 
necessary to perform boundary-layer solutions at the next T w entry of 1500° R 
for the current bracketing values of time and m* (Solutions 21 through 24) . 
The charring-ablation solution then recommenced and continued until the tabu- 
lar entry time of 4375 seconds, at which point a boundary-layer solution was 
obtained (Solution 25) . In order to proceed further with the charring-ablation 
solution, it was necessary to obtain the ablation temperatures (Solutions 26 

and 27) and current m* and T entries (Solutions 28 through 31) for the 

g w 

next time entry, 0^ = 4400 seconds . 

Returning to the transient charring-ablation solution, it can be seen 

that substantial pyrolysis was beginning to occur. When an m* of 0.1 was 
. , ^ 
attained, it was necessary to compute ablation temperatures (Solutions 32 and 

33) for the next m* entry of 0.2 for the two currently bracketing times, 

and to compute boundary layers for this new m* at the current T and times 

g w 

(Solutions 34 through 37) . Returning again to the charring-ablation solution, 
an m* of 0.2 was soon reached. Again, it was necessary to perform boundary- 
layer solutions to obtain ablation temperatures for the new m* of 0.3 (Solu- 
tions 38 and 39)* and to obtain boundary-layer solutions at the current 


Computer printout for boundary-layer Solution 39 is presented as Figure 3 of 
this report. 
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bracketing times and wall temperatures for this new m*. Boundary-layer 
Solution 40 was nonconvergent and the run was terminated at this point. 

The reason for this nonconvergent boundary-layer solution was that in- 
adequate convergence damping was applied for this particular problem. The 
combination of a high mass injection rate and low surface temperature pro- 
duced a large equilibrium concentration of CH 4 in the boundary-layer which, 
in turn, produced a pronounced temperature reversal, A convergent solution 
for this boundary-layer point has since been obtained by use of a different 
damping scheme. This convergence difficulty has been brought to light to 
demonstrate the fact that the boundary-layer procedure and in turn the chem- 
istry subroutines must be one hundred percent reliable for all potential 
flight conditions to be encountered in order to obtain a coupled solution 
without interruption. 

Temperature distributions through the charring— ablation material and 
boundary layer are presented in Table II for 9 = 4310, 4348 and 4375 seconds, 
respectively. Note that there is a substantial temperature change between 
the wall and the first node out into the boundary layer, suggesting that it 
might be advisable to add a few nodes or rearrange nodal spacing. 

The heat flux available to the charring-ablation program at an inter- 
mediate time is computed by linear interpolation of information derived from 
the boundary-layer solutions at the tabulated times. This is compared to the 
heat flux calculated by the boundary-layer procedure for the actual m*, m* 

and T in the following tabulation. It can be seen that the interpolated 
w 


Time, Sec 

Surface Heat Flux, Btu/sec ft 3 

From CMA Interpolation 
of BLIMP Calculations 

From Direct BLIMP Calculation 

4310 

0.154 

0.1541 

4348 

0.570 

0.5688 

4375 

1.58 

1.581 


values agree very closely with the actual values, indicating that the interpo- 
lation has been very accurate. 

The above-described coupled run, including 40 boundary-layer solutions, 
required 16 minutes and 33 seconds on the Univac 1108 computer. It is esti- 
mated that this is only 15 percent of the entire SA 202 reentry trajectory. 

At this rate a stagnation point solution for the entire trajectory would re- 
quire 1 hour and 50 minutes and the calculation of the nonsimilar boundary 
Taygjr around the Apollo reentry vehicle would require 18 hours and 20 minutes. 
This figure is misrepresentative for three reasons, two having to do with 
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TABLE II 


TEMPERATURE DISTRIBUTION FROM THE BACKWALL TO THE 
EDGE OF THE BOUNDARY LAYER FROM COUPLED SOLUTION: 
APOLLO STAGNATION POINT; SA 202 TRAJECTORY 


Nodal Point 


Nodal Temperatures, 
for Trajectory Times 

°R 

Of 

4,310 

sec 

4,348 sec 

4,375 sec 

Boundary Layer 





7 (edge) 

7,877 

8,584 

9, 228 

6 

7,504 

8,215 

8,858 

5 

6,995 

7,634 

8, 213 

4 

6,651 

7,226 

7,742 

3 

6,046 

6, 500 

6,891 

2 

3,493 

3,803 

4, 087 

1 (wall) 

530.0 

748.4 

1,099 

Charring Ablator 





1 (wall) 

530.0 

748.4 

1,099.0 

2 



737.4 

1, 070.3 

3 



717.2 

1,018.6 

4 



698.7 

971.8 

5 



681.7 

929.4 

6 



666.2 

890.8 

7 



651.9 

855.8 

8 



638.9 

824.1 

9 



627.1 

795.4 

10 



616.4 

769.4 

11 



606.7 

745.9 

12 



598.0 

724.7 

13 



590.1 

705.4 

14 



583.0 

688.0 

15 



576.7 

672.1 

16 



571.0 

657.8 

17 



564.7 

641.4 

18 



558.3 

624.6 

19 



553.0 

610.2 

20 



548.6 

597.8 

36 



530.2 

531.9 

37 

I 

530.1 

531.1 

38 (backwall) 

530.0 

530.0 

530.7 
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the number of chemistry iterations per boundary-layer iteration. In the pres- 
ent series* of 40 boundaty-layer solutions, an average of 14.3 iterations per 
solution was required. This is substantially greater than the 4 to 6 itera- 
tions usually required for a stagnation point solution, an uninspired uni- 
form damping factor of 0.6 having been applied to all corrections to improve 
reliability.* With an improved damping scheme it is reasonable to expect an 
average of 6 iterations per stagnation solution to be adequate. Since the 
boundary-layer calculations control the computational time, this would re- 
duce the time required to achieve a complete stagnation point history to 46 
minutes. Experience has shown that boundary-layer solutions for downstream 
stations converge very rapidly (about 2 to 4 iterations per solution) , and, 
furthermore, that the chemistry iteration is faster by a factor of two or so, 
both because the residual values for dependent variables provide good first 
guesses. This latter point is significant since the boundary-layer program 
spends about 70 percent of its time in the chemistry iteration. Applying 
these corrections to the original estimate# the computational time for the 
nine downstream stations could be expected to be reduced to 2 hours and 20 
minutes. Thus, the total computational time for a coupled transient non- 
similar solution around the Apollo heat-shield during an 800 second reentry 
trajectory would be approximately 3 hours on the Univac 1108. On the other 
side of the ledger, if one were to increase the number of nodal points across 
the boundary layer from 7 to 10 to insure accuracy# the computational time 
would increase to 4 hours and 30 minutes, since the computational time is 
approximately proportional to the number of nodal points with this number of 
nodes .** 

At the time of this writing, this is the only coupled transient charring 
ablation problem which has been attempted. It is significant, therefore, that 
the logic of the CABLE program (which calls the CMA and BLIMP programs as sub- 
routines) and the coupling mechanics are really quite straightforward. Further- 
more, the CMA program operates very nearly as it has on hundreds of occasions 
in the CMA/ ACE approach and the BLIMP program is also called upon to perform 
for the same classes of boundary conditions for which considerable experience 
has been gained. Therefore, it is fair to say that the CABLE program is rea- 
sonably well checked out. The major concerns are the improvement of computa- 
tional time and increase in reliability of the boundary layer and chemistry 
procedures. Recommendations for improvements are presented in Section 6. 


*This is very uneconomical in general but, as discussed previously, was in- 
adequate for the fortieth boundary-layer solution. 

**If a considerably larger number of nodes were to be used, say 15 or more, the 
computational time would become proportional to the cube of the number of 
nodal points, since matrix inversion would then be controlling. 
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SECTION 6 

SUMMARY, CONCLUSIONS AND RECOMMENDATIONS 

In the Parts II through VI of this series of reports and in this summary 
report (Part I) a substantial analytical effort has been described. In order 
to bring this effort into perspective, this section presents (1) a brief 
summary of the analytical studies and their product computer codes, (2) the 
conclusions which can be drawn relative to the methods adopted and their appli- 
caj °ility , and (3) recommendations for broader applications of the techniques 
and further analytical and numerical developments. 


6.1 SUMMARY 

The overall goal of this study has been the development of advanced theo- 
retical techniques and their implementation in the form of computer programs 
for the evaluation of ablation-material performance. The thorough treatment 
of the intimate thermal-chemical coupling between boundary- layer and ablative 
surface mechanisms has been the primary specific goal. 

The development of generalized and comprehensive predictive tools has 
been a dominant goal of the study, a factor which has necessitated the develop- 
ment of several novel analytical and numerical procedures. These procedures 
involve both the component parts of the analysis and the means whereby the 
parts are coupled. 

The general status of the procedures developed or extended during the 
current effort are summarized in Table I. The status of the analyses pre- 
sented in this series of reports is entered under "Model" in the table, whereas 
the status of the computer programs is indicated under "Operational Status." 

The table is directed specifically toward the Coupled Ablator/Boundary Layer/ 
Environment (CABLE) master program. This program utilizes three major compo- 
nents, namely, the evaluation of the in-depth response of the charring ablator 
and its backup materials (charring Material Ablation (CMA) program) , the 
characterization of the nonsimilar laminar multicomponent chemically-active 
boundary layer (Boundary Layer Integral Matrix (BLIMP) program) , and the 
evaluation of the chemical state in both open and closed systems (subroutine 
version of the Aerotherm Chemical Equilibrium (ACE) program) . The reader is 
referred to Table I for an overview of the developments. In the remainder of 
this section some of the highlights are reviewed. 

As a consequence of the ultimate coupling goal of this effort, computa- 
tional speed, reliability, compatibility and generality were significant fac- 
tors in the development of both theoretical and numerical techniques. Thus 
new boundary-layer and chemical-state procedures were developed. These have 
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demonstrated to be^major improvements (as judged by the above-listed factors) 
over existing procedures in all four areas. However, further improvements in 
computational speed and reliability of the boundary-layer and chemistry pro- 
grams are desirable and are discussed in Section 6.3 

The boundary layer is solved by an integral matrix solution procedure 
which yields accurate solutions with relatively few nodal points and thus 
relatively few entries into the conservation equations. Because of the major 
computational time associated with the evaluation of the state-dependent 
parameters in a general chemical system, speed is significantly enhanced by 
this minimization of nodal points. This numerical procedure possesses a high 
degree of mathematical formality and is applicable to a broad range of prob- 
lems with varied input and boundary conditions. An equilibrium chemistry 
model is incorporated into the computer programs which is applicable to all 
chemical systems. An analysis for extending this to general mixed equilibrium- 
nonequilibrium chemical systems is presented which should avoid the pitfalls 
often experienced as equilibrium is approached. Newly developed models for 
multicomponent transport properties including unequal diffusion and thermal 
diffusion coefficients for all species have demonstrated that accurate and 
time-saving approximations can be adopted with no significant loss of accu- 
racy. Finally, an implicit finite difference procedure is employed which 
accurately characterizes the in— depth response of the charring ablator. 

Two approaches have been developed for obtaining coupled solutions. The 
boundary layer and charring ablator numerical procedures can be coupled by 
use of the CABLE calling routine, or the CMA program can be operated indepen- 
dently utilizing input for the surface thermochemical boundary condition as 
provided by ACE (which represents the boundary layer by bulk transfer coef- 
ficients) . 

In the fully coupled procedure, the transient charring ablation response 
at a particular streamwise station is completed prior to advancing to the 
next station. Thus historic (i.e., upstream) boundary- layer information re- 
quired for a starting solution and for calculation of nonsimilar terms is 
stored, whereas no charring ablation information at stations other than the 
one under current consideration is stored. The approach which has been adopted 
has the primary advantages that (1) it provides an implicit solution while 
avoiding extrapolations of surface boundary conditions and iterative repeti- 
tion of charring ablation and boundary-layer solutions, (2) the CMA program 
operates very nearly the same as it does in the convective transfer coeffi- 
cient approach, the BLIMP program providing tabular input data as needed ( full 
tables being supplied a priori in the transfer coefficient approach) , and 
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(3) the BLIMP program operates with assigned wall conditions (a more direct 
option than the alternative of introducing overall wall mass and energy bal- 
ances into the boundary- layer solution) . 

6 . 2 CONCLUSIONS 

The following major conclusions can be drawn regarding the research de- 
scribed in this report: 

1- Theoretical analyses and computational techniques have been developed 
which are believed to extend substantially the state of the art for 
the characterization of: (a) nonsimilar, laminar, chemically- 

reacting boundary layers with unequal diffusion and thermal diffusion 
coefficients for all species and radiation absorption and emission; 

(b) mixed e qu i librium- none qui librium, homogeneous or heterogeneous 
chemical systems; and (c) the transient response of ablating mate- 
rials which decompose in depth, considering a general internal de- 
composition model (including coking kinetics) and with detailed 
consideration of chemical interactions at the surface. The specific 
models which have been employed are summarized in Table I. These 
procedures are self-consistent such as to permit full coupling. 

2. The computational techniques have been incorporated for the most 
part into a set of computer programs which can either be operated 
independently to yield solutions for the boundary layer, the chemical 
state, and/or the charring ablator, or be operated simultaneously 

to yield fully coupled solutions. The operational status of the 
program components are also summarized in Table I. 

3. These computer codes are applicable to all chemical systems and thus 
can be used to obtain solutions for any environment, ablation-material 
combination. These solutions should compare favorably with experi- 
mental data obtained in laminar flows as long as surface thermochemical 
considerations control the ablation process and if material and gas- 
phase properties are adequately known. 

4. For those material-environment combinations where subsurface chemis- 
try and/or surface mechanical removal mechanisms are not adequately 
described by the models currently employed, agreement with experi- 
mental data cannot be assured. 

6.3 RECOMMENDATIONS 

Recommendations for further work can be divided into three general cate- 
gories; recommendations on the use of the computer programs in their present 
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form, recommendations for improving certain computational aspects of the solu- 
tion procedures , and recommendations pertaining to the inclusion of additional 
physicochemical phenomena in the computer codes . These three groups are dis- 
cussed in the following subsections. 

6.3,1 Program Employment 

As mentioned previously, the programs in their present form are adequate 
for many materials and environments, in particular those where surface thermo- 
chemical considerations control the ablation process and where material and 
gas-phase properties are adequately known. When this is the case, utilization 
of the programs is straightforward. However, when other considerations are 
significant, some "education" of the programs may be required for the specific 
material (s) of interest. The following approach is recommended to accomplish 
this objective: 

1. Perform transient charring ablation calculations with the CMA program. 
Option 2 (assigned surface temperatures and ablation rates as measured 
during static test operation) and correlate with experimental internal 
temperature distributions as a check on thermal property data. Post- 
test chemical and physical analyses, in-depth, might also indicate 
the necessity for inclusion of such mechanisms as coking or erosion 
within the char layer or condensed phase reactions in-depth. 

2. Perform transient charring ablation calculations with CMA/EST or 
CMA/ACE programs (which represent the boundary layer by bulk trans- 
fer coefficients) after first obtaining the necessary transfer coef- 
ficients from calorimetric data, from separate operation of the 
BLIMP program, or if possible, from cross-correlation of the results 
of both methods. With these calculations, attempt to correlate sur- 
face recession and surface temperature with experimental data. To 
achieve satisfactory correlation it may be necessary to extend these 
programs to include such considerations as the allowance of multiple 
condensed phase surface materials (e.g., liquid droplets on an other- 
wise carbon surface, or liquid solutions) , and surface reaction 
kinetics. The latter can be taken into account by construction of a 
specialized KINET subroutine for the ACE program. The proper intro- 
duction of any of the mechanisms mentioned above is dependent on both 
the formulation of valid models and the experimental or theoretical 
evaluation of the appropriate properties. 

3. Once reasonably good correlation is obtained, incorporate any changes 
in the ACE (or ACE-KINET) program into the surface chemistry routines 
of the BLIMP subroutine to the CABLE program, and utilize the CABLE 
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program to obtain coupled nonsimilar solutions. Again, the predic- 
tions should be correlated with experimental data. 

In the event that a large number of calculations are to be made, paramet- 
ric solutions using the CABLE program should be generated. The results should 
then be used to develop formulations for convective transfer coefficients spe- 
cific to the particular materials and applications of interest for use with 
the more economical and more expedient CMA/ACE or CMA/ACE-KINET programs. 

6.3.2 Improved Computational Details 

The CMA program for in-depth response does not appear to require any ex- 
tensive upgrading in computational details. For rapid transients, during 
which it is desirable to reduce the time-step size (in order to follow the 
transient more accurately) to values below the maximum which would be allowed 
by an explicit solution, it would be of value to be able to switch to a full- 
explicit in-depth energy scheme in order to minimize computation time. This 
Could be simply done. For transients involving extremely high heat fluxes, 
on the other hand, it appears that neither a full-explicit scheme nor the 
mixed implicit-energy, explicit-density scheme presently used in the CMA pro- 
gram is always adequate. For these transients, both the energy equation and 
the decomposition equation appear to require implicit treatment to prevent 
the solution from disintegrating. This problem deserves further attention. 

With regard to the BLIMP program, it is recommended that effort be di- 
rected to the development of a general convergence damping scheme. In those 
instances where inexact analytic derivatives are currently employed in the 
Newton-Raphson iteration process, it may be desirable to compute exact analyt- 
ical derivatives or to approximate at least some of these derivatives by 
finite-difference relations. These are especially important if the BLIMP 
program is to be used extensively as a subroutine to the CABLE program be- 
cause 100 percent reliability is required. Finally, for the purpose of com- 
putational economy, a streamlined CABLE program should be developed where 
considerably fewer boundary-layer solutions are required, interpolation be- 
tween these solutions being guided by the considerably faster film-coefficient 
approach . 

6.3.3 Additions to the Physical Model 

Since the CMA in-depth computation includes only the physics of the basic 
pyrolysis problem, specific materials with important additional subsurface 
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events are not all well described by the program. Physical events here include 
coking and subsurface char erosion due to interaction with the pyrolysis gas, 
chemical kinetics of pyrolysis-gas cracking reactions as the gas flows through 
the char, thick liquid layer run-off (the present chemical programs account 
for only thin, nonviscous liquid-layer removal) , additional subsurface chemical 
reactions such as carbon-silica reactions, thermal expansion effects, and me- 
chanical damage to weak chars. 

For many materials none of these additional effects is of interest. For 
other materials, the relative importance of these additional events has gen- 
erally not been determined. For such materials the probable importance of the 
various physical events should be defined through well-planned experimentation 
and post-test model examination coupled with performance prediction with the 
existing computer programs.. Such an approach has been outlined in Section 
6.3.1. Only after such experimentation and study can additional analysis and 
programming efforts be undertaken with real confidence. In this regard, an 
approach for considering coking reactions in-depth, presented in Part VI of 
this series of reports, will be useful in the event that char density buildup 
appears to be an important process. A rather extensive programming effort 
would be required to incorporate this model into the CMA program. 

The assessment of the applicability of the boundary-layer procedure in 
its present operational form for a particular ablation material in a specific 
application is much more straightforward since it is possible to estimate 
under what conditions phenomena not currently considered come into play. The 
following computer program developments may be required depending primarily 
upon the flight conditions of interest: extension to a turbulent boundary 

layer, to include radiation absorption and emission within the boundary layer, 
to include general nonequilibrium, to include an entropy layer, or to allow 
a nonadiabatic inviscid flow field. The effort required to accomplish each 
of these extensions is discussed briefly in the remaining paragraphs. 

The boundary layer computational procedure could be extended to turbulent 
boundary layers by the use of eddy transport properties based, for example, on 
the laws of the wall and the wake adopted in Ref. 21. Radiation absorp- 
tion and emission could be included by the use of the one-dimensional model 
presented in Section 2.1.3 and discussed in more detail in Appendix E to 
Part III. 

The first step to the inclusion of nonequilibrium effects within the 
boundary layer could be to perform boundary-layer calculations while consider- 
ing the kinetically controlled reactions to be frozen within the boundary layer, 
but including these reactions at the boundary-layer edge and including surface 
catalyzed and heterogeneous reactions at the wall. All equilibrium reactions 
would then be included during all aspects of the calculations. Subsequently 
nonequilibrium streamtube calculations would then be performed to obtain the 
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complete nonequilibrium flow field. This is permissible as a first approxi- 
mation since ablation rates are generally fairly insensitive to nonequilibrium 
effects within the boundary layer. The extension of the boundary-layer pro- 
cedure to include the general mixed equilibrium-nonequilibrium model presented 
in Section 4 and described in more detail in Part V of this series of reports 
is believed to be practical but would require a rather extensive programming 
effort . 

The BLIMP program is presently programmed to include entropy-layer effects 
but this option has never been activated. Consideration of a nonadiabatic flow 
field would require the reorganization of some program logic, but this could 
be done without major difficulty. 
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